File indexing completed on 2022-02-21 23:14:04
0001 #include "FWCore/Framework/interface/Frameworkfwd.h"
0002 #include "FWCore/Framework/interface/stream/EDProducer.h"
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "FWCore/Framework/interface/EventSetup.h"
0005 #include "FWCore/Framework/interface/ESHandle.h"
0006 #include "FWCore/Framework/interface/ConsumesCollector.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008
0009 #include "RecoMTD/DetLayers/interface/MTDDetLayerGeometry.h"
0010 #include "RecoMTD/Records/interface/MTDRecoGeometryRecord.h"
0011
0012 #include "MagneticField/Engine/interface/MagneticField.h"
0013 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0014
0015 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0016
0017 #include "DataFormats/TrackerRecHit2D/interface/MTDTrackingRecHit.h"
0018
0019 #include "RecoMTD/DetLayers/interface/MTDTrayBarrelLayer.h"
0020 #include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"
0021
0022 #include "DataFormats/ForwardDetId/interface/BTLDetId.h"
0023 #include "DataFormats/ForwardDetId/interface/ETLDetId.h"
0024 #include "DataFormats/ForwardDetId/interface/MTDChannelIdentifier.h"
0025 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0026
0027 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0028 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
0029
0030 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0031 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0032 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0033
0034 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0035
0036 #include "RecoMTD/TransientTrackingRecHit/interface/MTDTransientTrackingRecHitBuilder.h"
0037 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0038
0039 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0040 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0041
0042 #include "TrackingTools/PatternTools/interface/TSCBLBuilderWithPropagator.h"
0043
0044 #include "RecoTracker/TransientTrackingRecHit/interface/Traj2TrackHits.h"
0045 #include "TrackingTools/TrackRefitter/interface/TrackTransformer.h"
0046
0047 #include <sstream>
0048
0049 #include "Geometry/CommonTopologies/interface/Topology.h"
0050
0051 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0052 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0053
0054 #include "DataFormats/Math/interface/GeantUnits.h"
0055 #include "DataFormats/Math/interface/LorentzVector.h"
0056 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0057 #include "DataFormats/Math/interface/Rounding.h"
0058
0059 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0060 #include "DataFormats/VertexReco/interface/Vertex.h"
0061
0062 using namespace std;
0063 using namespace edm;
0064 using namespace reco;
0065
0066 namespace {
0067 constexpr float c_cm_ns = geant_units::operators::convertMmToCm(CLHEP::c_light);
0068 constexpr float c_inv = 1.0f / c_cm_ns;
0069
0070 class MTDHitMatchingInfo {
0071 public:
0072 MTDHitMatchingInfo() {
0073 hit = nullptr;
0074 estChi2 = std::numeric_limits<float>::max();
0075 timeChi2 = std::numeric_limits<float>::max();
0076 }
0077
0078
0079 inline bool operator<(const MTDHitMatchingInfo& m2) const {
0080
0081 constexpr float chi2_cut = 10.f;
0082 constexpr float low_weight = 3.f;
0083 constexpr float high_weight = 8.f;
0084 if (timeChi2 < chi2_cut && m2.timeChi2 < chi2_cut)
0085 return chi2(low_weight) < m2.chi2(low_weight);
0086 else
0087 return chi2(high_weight) < m2.chi2(high_weight);
0088 }
0089
0090 inline float chi2(float timeWeight = 1.f) const { return estChi2 + timeWeight * timeChi2; }
0091
0092 const MTDTrackingRecHit* hit;
0093 float estChi2;
0094 float timeChi2;
0095 };
0096
0097 class TrackSegments {
0098 public:
0099 TrackSegments() = default;
0100
0101 inline uint32_t addSegment(float tPath, float tMom2) {
0102 segmentPathOvc_.emplace_back(tPath * c_inv);
0103 segmentMom2_.emplace_back(tMom2);
0104 nSegment_++;
0105
0106 LogTrace("TrackExtenderWithMTD") << "addSegment # " << nSegment_ << " s = " << tPath
0107 << " p = " << std::sqrt(tMom2);
0108
0109 return nSegment_;
0110 }
0111
0112 inline float computeTof(float mass_inv2) const {
0113 float tof(0.f);
0114 for (uint32_t iSeg = 0; iSeg < nSegment_; iSeg++) {
0115 float gammasq = 1.f + segmentMom2_[iSeg] * mass_inv2;
0116 float beta = std::sqrt(1.f - 1.f / gammasq);
0117 tof += segmentPathOvc_[iSeg] / beta;
0118
0119 LogTrace("TrackExtenderWithMTD") << " TOF Segment # " << iSeg + 1 << " p = " << std::sqrt(segmentMom2_[iSeg])
0120 << " tof = " << tof;
0121 }
0122
0123 return tof;
0124 }
0125
0126 inline uint32_t size() const { return nSegment_; }
0127
0128 inline uint32_t removeFirstSegment() {
0129 if (nSegment_ > 0) {
0130 segmentPathOvc_.erase(segmentPathOvc_.begin());
0131 segmentMom2_.erase(segmentMom2_.begin());
0132 nSegment_--;
0133 }
0134 return nSegment_;
0135 }
0136
0137 inline std::pair<float, float> getSegmentPathAndMom2(uint32_t iSegment) const {
0138 if (iSegment >= nSegment_) {
0139 throw cms::Exception("TrackExtenderWithMTD") << "Requesting non existing track segment #" << iSegment;
0140 }
0141 return std::make_pair(segmentPathOvc_[iSegment], segmentMom2_[iSegment]);
0142 }
0143
0144 uint32_t nSegment_ = 0;
0145 std::vector<float> segmentPathOvc_;
0146 std::vector<float> segmentMom2_;
0147 };
0148
0149 struct TrackTofPidInfo {
0150 float tmtd;
0151 float tmtderror;
0152 float pathlength;
0153
0154 float betaerror;
0155
0156 float dt;
0157 float dterror;
0158 float dtchi2;
0159
0160 float dt_best;
0161 float dterror_best;
0162 float dtchi2_best;
0163
0164 float gammasq_pi;
0165 float beta_pi;
0166 float dt_pi;
0167
0168 float gammasq_k;
0169 float beta_k;
0170 float dt_k;
0171
0172 float gammasq_p;
0173 float beta_p;
0174 float dt_p;
0175
0176 float prob_pi;
0177 float prob_k;
0178 float prob_p;
0179 };
0180
0181 enum class TofCalc { kCost = 1, kSegm = 2, kMixd = 3 };
0182
0183 const TrackTofPidInfo computeTrackTofPidInfo(float magp2,
0184 float length,
0185 TrackSegments trs,
0186 float t_mtd,
0187 float t_mtderr,
0188 float t_vtx,
0189 float t_vtx_err,
0190 bool addPIDError = true,
0191 TofCalc choice = TofCalc::kCost) {
0192 constexpr float m_pi = 0.13957018f;
0193 constexpr float m_pi_inv2 = 1.0f / m_pi / m_pi;
0194 constexpr float m_k = 0.493677f;
0195 constexpr float m_k_inv2 = 1.0f / m_k / m_k;
0196 constexpr float m_p = 0.9382720813f;
0197 constexpr float m_p_inv2 = 1.0f / m_p / m_p;
0198
0199 TrackTofPidInfo tofpid;
0200
0201 tofpid.tmtd = t_mtd;
0202 tofpid.tmtderror = t_mtderr;
0203 tofpid.pathlength = length;
0204
0205 auto deltat = [&](const float mass_inv2, const float betatmp) {
0206 float res(1.f);
0207 switch (choice) {
0208 case TofCalc::kCost:
0209 res = tofpid.pathlength / betatmp * c_inv;
0210 break;
0211 case TofCalc::kSegm:
0212 res = trs.computeTof(mass_inv2);
0213 break;
0214 case TofCalc::kMixd:
0215 res = trs.computeTof(mass_inv2) + tofpid.pathlength / betatmp * c_inv;
0216 break;
0217 }
0218 return res;
0219 };
0220
0221 tofpid.gammasq_pi = 1.f + magp2 * m_pi_inv2;
0222 tofpid.beta_pi = std::sqrt(1.f - 1.f / tofpid.gammasq_pi);
0223 tofpid.dt_pi = deltat(m_pi_inv2, tofpid.beta_pi);
0224
0225 tofpid.gammasq_k = 1.f + magp2 * m_k_inv2;
0226 tofpid.beta_k = std::sqrt(1.f - 1.f / tofpid.gammasq_k);
0227 tofpid.dt_k = deltat(m_k_inv2, tofpid.beta_k);
0228
0229 tofpid.gammasq_p = 1.f + magp2 * m_p_inv2;
0230 tofpid.beta_p = std::sqrt(1.f - 1.f / tofpid.gammasq_p);
0231 tofpid.dt_p = deltat(m_p_inv2, tofpid.beta_p);
0232
0233 tofpid.dt = tofpid.tmtd - tofpid.dt_pi - t_vtx;
0234 tofpid.dterror = sqrt(tofpid.tmtderror * tofpid.tmtderror + t_vtx_err * t_vtx_err);
0235 tofpid.betaerror = 0.f;
0236 if (addPIDError) {
0237 tofpid.dterror =
0238 sqrt(tofpid.dterror * tofpid.dterror + (tofpid.dt_p - tofpid.dt_pi) * (tofpid.dt_p - tofpid.dt_pi));
0239 tofpid.betaerror = tofpid.beta_p - tofpid.beta_pi;
0240 }
0241
0242 tofpid.dtchi2 = (tofpid.dt * tofpid.dt) / (tofpid.dterror * tofpid.dterror);
0243
0244 tofpid.dt_best = tofpid.dt;
0245 tofpid.dterror_best = tofpid.dterror;
0246 tofpid.dtchi2_best = tofpid.dtchi2;
0247
0248 tofpid.prob_pi = -1.f;
0249 tofpid.prob_k = -1.f;
0250 tofpid.prob_p = -1.f;
0251
0252 if (!addPIDError) {
0253
0254 float chi2_pi = tofpid.dtchi2;
0255 float chi2_k =
0256 (tofpid.tmtd - tofpid.dt_k - t_vtx) * (tofpid.tmtd - tofpid.dt_k - t_vtx) / (tofpid.dterror * tofpid.dterror);
0257 float chi2_p =
0258 (tofpid.tmtd - tofpid.dt_p - t_vtx) * (tofpid.tmtd - tofpid.dt_p - t_vtx) / (tofpid.dterror * tofpid.dterror);
0259
0260 float rawprob_pi = exp(-0.5f * chi2_pi);
0261 float rawprob_k = exp(-0.5f * chi2_k);
0262 float rawprob_p = exp(-0.5f * chi2_p);
0263 float normprob = 1.f / (rawprob_pi + rawprob_k + rawprob_p);
0264
0265 tofpid.prob_pi = rawprob_pi * normprob;
0266 tofpid.prob_k = rawprob_k * normprob;
0267 tofpid.prob_p = rawprob_p * normprob;
0268
0269 float prob_heavy = 1.f - tofpid.prob_pi;
0270 constexpr float heavy_threshold = 0.75f;
0271
0272 if (prob_heavy > heavy_threshold) {
0273 if (chi2_k < chi2_p) {
0274 tofpid.dt_best = (tofpid.tmtd - tofpid.dt_k - t_vtx);
0275 tofpid.dtchi2_best = chi2_k;
0276 } else {
0277 tofpid.dt_best = (tofpid.tmtd - tofpid.dt_p - t_vtx);
0278 tofpid.dtchi2_best = chi2_p;
0279 }
0280 }
0281 }
0282 return tofpid;
0283 }
0284
0285 bool getTrajectoryStateClosestToBeamLine(const Trajectory& traj,
0286 const reco::BeamSpot& bs,
0287 const Propagator* thePropagator,
0288 TrajectoryStateClosestToBeamLine& tscbl) {
0289
0290 TrajectoryStateOnSurface stateForProjectionToBeamLineOnSurface =
0291 traj.closestMeasurement(GlobalPoint(bs.x0(), bs.y0(), bs.z0())).updatedState();
0292
0293 if (!stateForProjectionToBeamLineOnSurface.isValid()) {
0294 edm::LogError("CannotPropagateToBeamLine") << "the state on the closest measurement isnot valid. skipping track.";
0295 return false;
0296 }
0297
0298 const FreeTrajectoryState& stateForProjectionToBeamLine = *stateForProjectionToBeamLineOnSurface.freeState();
0299
0300 TSCBLBuilderWithPropagator tscblBuilder(*thePropagator);
0301 tscbl = tscblBuilder(stateForProjectionToBeamLine, bs);
0302
0303 return tscbl.isValid();
0304 }
0305
0306 bool trackPathLength(const Trajectory& traj,
0307 const TrajectoryStateClosestToBeamLine& tscbl,
0308 const Propagator* thePropagator,
0309 float& pathlength,
0310 TrackSegments& trs) {
0311 pathlength = 0.f;
0312
0313 bool validpropagation = true;
0314 float oldp = traj.measurements().begin()->updatedState().globalMomentum().mag();
0315 float pathlength1 = 0.f;
0316 float pathlength2 = 0.f;
0317
0318
0319 for (auto it = traj.measurements().begin(); it != traj.measurements().end() - 1; ++it) {
0320 const auto& propresult = thePropagator->propagateWithPath(it->updatedState(), (it + 1)->updatedState().surface());
0321 float layerpathlength = std::abs(propresult.second);
0322 if (layerpathlength == 0.f) {
0323 validpropagation = false;
0324 }
0325 pathlength1 += layerpathlength;
0326 trs.addSegment(layerpathlength, (it + 1)->updatedState().globalMomentum().mag2());
0327 LogTrace("TrackExtenderWithMTD") << "TSOS " << std::fixed << std::setw(4) << trs.size() << " R_i " << std::fixed
0328 << std::setw(14) << it->updatedState().globalPosition().perp() << " z_i "
0329 << std::fixed << std::setw(14) << it->updatedState().globalPosition().z()
0330 << " R_e " << std::fixed << std::setw(14)
0331 << (it + 1)->updatedState().globalPosition().perp() << " z_e " << std::fixed
0332 << std::setw(14) << (it + 1)->updatedState().globalPosition().z() << " p "
0333 << std::fixed << std::setw(14) << (it + 1)->updatedState().globalMomentum().mag()
0334 << " dp " << std::fixed << std::setw(14)
0335 << (it + 1)->updatedState().globalMomentum().mag() - oldp;
0336 oldp = (it + 1)->updatedState().globalMomentum().mag();
0337 }
0338
0339
0340 auto const& tscblPCA = tscbl.trackStateAtPCA();
0341 auto const& aSurface = traj.direction() == alongMomentum ? traj.firstMeasurement().updatedState().surface()
0342 : traj.lastMeasurement().updatedState().surface();
0343 pathlength2 = thePropagator->propagateWithPath(tscblPCA, aSurface).second;
0344 if (pathlength2 == 0.f) {
0345 validpropagation = false;
0346 }
0347 pathlength = pathlength1 + pathlength2;
0348 trs.addSegment(pathlength2, tscblPCA.momentum().mag2());
0349 LogTrace("TrackExtenderWithMTD") << "TSOS " << std::fixed << std::setw(4) << trs.size() << " R_e " << std::fixed
0350 << std::setw(14) << tscblPCA.position().perp() << " z_e " << std::fixed
0351 << std::setw(14) << tscblPCA.position().z() << " p " << std::fixed << std::setw(14)
0352 << tscblPCA.momentum().mag() << " dp " << std::fixed << std::setw(14)
0353 << tscblPCA.momentum().mag() - oldp;
0354 return validpropagation;
0355 }
0356
0357 bool trackPathLength(const Trajectory& traj,
0358 const reco::BeamSpot& bs,
0359 const Propagator* thePropagator,
0360 float& pathlength,
0361 TrackSegments& trs) {
0362 pathlength = 0.f;
0363
0364 TrajectoryStateClosestToBeamLine tscbl;
0365 bool tscbl_status = getTrajectoryStateClosestToBeamLine(traj, bs, thePropagator, tscbl);
0366
0367 if (!tscbl_status)
0368 return false;
0369
0370 return trackPathLength(traj, tscbl, thePropagator, pathlength, trs);
0371 }
0372
0373 }
0374
0375 template <class TrackCollection>
0376 class TrackExtenderWithMTDT : public edm::stream::EDProducer<> {
0377 public:
0378 typedef typename TrackCollection::value_type TrackType;
0379 typedef edm::View<TrackType> InputCollection;
0380
0381 TrackExtenderWithMTDT(const ParameterSet& pset);
0382
0383 template <class H, class T>
0384 void fillValueMap(edm::Event& iEvent, const H& handle, const std::vector<T>& vec, const edm::EDPutToken& token) const;
0385
0386 void produce(edm::Event& ev, const edm::EventSetup& es) final;
0387
0388 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0389
0390 TransientTrackingRecHit::ConstRecHitContainer tryBTLLayers(const TrajectoryStateOnSurface&,
0391 const Trajectory& traj,
0392 const float,
0393 const float,
0394 const TrackSegments&,
0395 const MTDTrackingDetSetVector&,
0396 const MTDDetLayerGeometry*,
0397 const MagneticField* field,
0398 const Propagator* prop,
0399 const reco::BeamSpot& bs,
0400 const float vtxTime,
0401 const bool matchVertex,
0402 MTDHitMatchingInfo& bestHit) const;
0403
0404 TransientTrackingRecHit::ConstRecHitContainer tryETLLayers(const TrajectoryStateOnSurface&,
0405 const Trajectory& traj,
0406 const float,
0407 const float,
0408 const TrackSegments&,
0409 const MTDTrackingDetSetVector&,
0410 const MTDDetLayerGeometry*,
0411 const MagneticField* field,
0412 const Propagator* prop,
0413 const reco::BeamSpot& bs,
0414 const float vtxTime,
0415 const bool matchVertex,
0416 MTDHitMatchingInfo& bestHit) const;
0417
0418 void fillMatchingHits(const DetLayer*,
0419 const TrajectoryStateOnSurface&,
0420 const Trajectory&,
0421 const float,
0422 const float,
0423 const TrackSegments&,
0424 const MTDTrackingDetSetVector&,
0425 const Propagator*,
0426 const reco::BeamSpot&,
0427 const float&,
0428 const bool,
0429 TransientTrackingRecHit::ConstRecHitContainer&,
0430 MTDHitMatchingInfo&) const;
0431
0432 RefitDirection::GeometricalDirection checkRecHitsOrdering(
0433 TransientTrackingRecHit::ConstRecHitContainer const& recHits) const {
0434 if (!recHits.empty()) {
0435 GlobalPoint first = gtg_->idToDet(recHits.front()->geographicalId())->position();
0436 GlobalPoint last = gtg_->idToDet(recHits.back()->geographicalId())->position();
0437
0438
0439 auto rFirst = first.mag2();
0440 auto rLast = last.mag2();
0441 if (rFirst < rLast)
0442 return RefitDirection::insideOut;
0443 if (rFirst > rLast)
0444 return RefitDirection::outsideIn;
0445 }
0446 LogDebug("TrackExtenderWithMTD") << "Impossible to determine the rechits order" << endl;
0447 return RefitDirection::undetermined;
0448 }
0449
0450 reco::Track buildTrack(const reco::TrackRef&,
0451 const Trajectory&,
0452 const Trajectory&,
0453 const reco::BeamSpot&,
0454 const MagneticField* field,
0455 const Propagator* prop,
0456 bool hasMTD,
0457 float& pathLength,
0458 float& tmtdOut,
0459 float& sigmatmtdOut,
0460 float& tofpi,
0461 float& tofk,
0462 float& tofp) const;
0463 reco::TrackExtra buildTrackExtra(const Trajectory& trajectory) const;
0464
0465 string dumpLayer(const DetLayer* layer) const;
0466
0467 private:
0468 edm::EDPutToken btlMatchChi2Token_;
0469 edm::EDPutToken etlMatchChi2Token_;
0470 edm::EDPutToken btlMatchTimeChi2Token_;
0471 edm::EDPutToken etlMatchTimeChi2Token_;
0472 edm::EDPutToken npixBarrelToken_;
0473 edm::EDPutToken npixEndcapToken_;
0474 edm::EDPutToken pOrigTrkToken_;
0475 edm::EDPutToken betaOrigTrkToken_;
0476 edm::EDPutToken t0OrigTrkToken_;
0477 edm::EDPutToken sigmat0OrigTrkToken_;
0478 edm::EDPutToken pathLengthOrigTrkToken_;
0479 edm::EDPutToken tmtdOrigTrkToken_;
0480 edm::EDPutToken sigmatmtdOrigTrkToken_;
0481 edm::EDPutToken tofpiOrigTrkToken_;
0482 edm::EDPutToken tofkOrigTrkToken_;
0483 edm::EDPutToken tofpOrigTrkToken_;
0484 edm::EDPutToken assocOrigTrkToken_;
0485
0486 edm::EDGetTokenT<InputCollection> tracksToken_;
0487 edm::EDGetTokenT<TrajTrackAssociationCollection> trajTrackAToken_;
0488 edm::EDGetTokenT<MTDTrackingDetSetVector> hitsToken_;
0489 edm::EDGetTokenT<reco::BeamSpot> bsToken_;
0490 edm::EDGetTokenT<GlobalPoint> genVtxPositionToken_;
0491 edm::EDGetTokenT<float> genVtxTimeToken_;
0492 edm::EDGetTokenT<VertexCollection> vtxToken_;
0493
0494 const bool updateTraj_, updateExtra_, updatePattern_;
0495 const std::string mtdRecHitBuilder_, propagator_, transientTrackBuilder_;
0496 std::unique_ptr<MeasurementEstimator> theEstimator;
0497 std::unique_ptr<TrackTransformer> theTransformer;
0498 edm::ESHandle<TransientTrackBuilder> builder_;
0499 edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> builderToken_;
0500 edm::ESHandle<TransientTrackingRecHitBuilder> hitbuilder_;
0501 edm::ESGetToken<TransientTrackingRecHitBuilder, TransientRecHitRecord> hitbuilderToken_;
0502 edm::ESHandle<GlobalTrackingGeometry> gtg_;
0503 edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> gtgToken_;
0504
0505 edm::ESGetToken<MTDDetLayerGeometry, MTDRecoGeometryRecord> dlgeoToken_;
0506 edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magfldToken_;
0507 edm::ESGetToken<Propagator, TrackingComponentsRecord> propToken_;
0508 edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> ttopoToken_;
0509
0510 const float estMaxChi2_;
0511 const float estMaxNSigma_;
0512 const float btlChi2Cut_;
0513 const float btlTimeChi2Cut_;
0514 const float etlChi2Cut_;
0515 const float etlTimeChi2Cut_;
0516
0517 const bool useVertex_;
0518 const bool useSimVertex_;
0519 const float dzCut_;
0520 const float bsTimeSpread_;
0521 };
0522
0523 template <class TrackCollection>
0524 TrackExtenderWithMTDT<TrackCollection>::TrackExtenderWithMTDT(const ParameterSet& iConfig)
0525 : tracksToken_(consumes<InputCollection>(iConfig.getParameter<edm::InputTag>("tracksSrc"))),
0526 trajTrackAToken_(consumes<TrajTrackAssociationCollection>(iConfig.getParameter<edm::InputTag>("trjtrkAssSrc"))),
0527 hitsToken_(consumes<MTDTrackingDetSetVector>(iConfig.getParameter<edm::InputTag>("hitsSrc"))),
0528 bsToken_(consumes<reco::BeamSpot>(iConfig.getParameter<edm::InputTag>("beamSpotSrc"))),
0529 updateTraj_(iConfig.getParameter<bool>("updateTrackTrajectory")),
0530 updateExtra_(iConfig.getParameter<bool>("updateTrackExtra")),
0531 updatePattern_(iConfig.getParameter<bool>("updateTrackHitPattern")),
0532 mtdRecHitBuilder_(iConfig.getParameter<std::string>("MTDRecHitBuilder")),
0533 propagator_(iConfig.getParameter<std::string>("Propagator")),
0534 transientTrackBuilder_(iConfig.getParameter<std::string>("TransientTrackBuilder")),
0535 estMaxChi2_(iConfig.getParameter<double>("estimatorMaxChi2")),
0536 estMaxNSigma_(iConfig.getParameter<double>("estimatorMaxNSigma")),
0537 btlChi2Cut_(iConfig.getParameter<double>("btlChi2Cut")),
0538 btlTimeChi2Cut_(iConfig.getParameter<double>("btlTimeChi2Cut")),
0539 etlChi2Cut_(iConfig.getParameter<double>("etlChi2Cut")),
0540 etlTimeChi2Cut_(iConfig.getParameter<double>("etlTimeChi2Cut")),
0541 useVertex_(iConfig.getParameter<bool>("useVertex")),
0542 useSimVertex_(iConfig.getParameter<bool>("useSimVertex")),
0543 dzCut_(iConfig.getParameter<double>("dZCut")),
0544 bsTimeSpread_(iConfig.getParameter<double>("bsTimeSpread")) {
0545 if (useVertex_) {
0546 if (useSimVertex_) {
0547 genVtxPositionToken_ = consumes<GlobalPoint>(iConfig.getParameter<edm::InputTag>("genVtxPositionSrc"));
0548 genVtxTimeToken_ = consumes<float>(iConfig.getParameter<edm::InputTag>("genVtxTimeSrc"));
0549 } else
0550 vtxToken_ = consumes<VertexCollection>(iConfig.getParameter<edm::InputTag>("vtxSrc"));
0551 }
0552
0553 theEstimator = std::make_unique<Chi2MeasurementEstimator>(estMaxChi2_, estMaxNSigma_);
0554 theTransformer = std::make_unique<TrackTransformer>(iConfig.getParameterSet("TrackTransformer"), consumesCollector());
0555
0556 btlMatchChi2Token_ = produces<edm::ValueMap<float>>("btlMatchChi2");
0557 etlMatchChi2Token_ = produces<edm::ValueMap<float>>("etlMatchChi2");
0558 btlMatchTimeChi2Token_ = produces<edm::ValueMap<float>>("btlMatchTimeChi2");
0559 etlMatchTimeChi2Token_ = produces<edm::ValueMap<float>>("etlMatchTimeChi2");
0560 npixBarrelToken_ = produces<edm::ValueMap<int>>("npixBarrel");
0561 npixEndcapToken_ = produces<edm::ValueMap<int>>("npixEndcap");
0562 pOrigTrkToken_ = produces<edm::ValueMap<float>>("generalTrackp");
0563 betaOrigTrkToken_ = produces<edm::ValueMap<float>>("generalTrackBeta");
0564 t0OrigTrkToken_ = produces<edm::ValueMap<float>>("generalTrackt0");
0565 sigmat0OrigTrkToken_ = produces<edm::ValueMap<float>>("generalTracksigmat0");
0566 pathLengthOrigTrkToken_ = produces<edm::ValueMap<float>>("generalTrackPathLength");
0567 tmtdOrigTrkToken_ = produces<edm::ValueMap<float>>("generalTracktmtd");
0568 sigmatmtdOrigTrkToken_ = produces<edm::ValueMap<float>>("generalTracksigmatmtd");
0569 tofpiOrigTrkToken_ = produces<edm::ValueMap<float>>("generalTrackTofPi");
0570 tofkOrigTrkToken_ = produces<edm::ValueMap<float>>("generalTrackTofK");
0571 tofpOrigTrkToken_ = produces<edm::ValueMap<float>>("generalTrackTofP");
0572 assocOrigTrkToken_ = produces<edm::ValueMap<int>>("generalTrackassoc");
0573
0574 builderToken_ = esConsumes<TransientTrackBuilder, TransientTrackRecord>(edm::ESInputTag("", transientTrackBuilder_));
0575 hitbuilderToken_ =
0576 esConsumes<TransientTrackingRecHitBuilder, TransientRecHitRecord>(edm::ESInputTag("", mtdRecHitBuilder_));
0577 gtgToken_ = esConsumes<GlobalTrackingGeometry, GlobalTrackingGeometryRecord>();
0578 dlgeoToken_ = esConsumes<MTDDetLayerGeometry, MTDRecoGeometryRecord>();
0579 magfldToken_ = esConsumes<MagneticField, IdealMagneticFieldRecord>();
0580 propToken_ = esConsumes<Propagator, TrackingComponentsRecord>(edm::ESInputTag("", propagator_));
0581 ttopoToken_ = esConsumes<TrackerTopology, TrackerTopologyRcd>();
0582
0583 produces<edm::OwnVector<TrackingRecHit>>();
0584 produces<reco::TrackExtraCollection>();
0585 produces<TrackCollection>();
0586 }
0587
0588 template <class TrackCollection>
0589 void TrackExtenderWithMTDT<TrackCollection>::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0590 edm::ParameterSetDescription desc, transDesc;
0591 desc.add<edm::InputTag>("tracksSrc", edm::InputTag("generalTracks"));
0592 desc.add<edm::InputTag>("trjtrkAssSrc", edm::InputTag("generalTracks"));
0593 desc.add<edm::InputTag>("hitsSrc", edm::InputTag("mtdTrackingRecHits"));
0594 desc.add<edm::InputTag>("beamSpotSrc", edm::InputTag("offlineBeamSpot"));
0595 desc.add<edm::InputTag>("genVtxPositionSrc", edm::InputTag("genParticles:xyz0"));
0596 desc.add<edm::InputTag>("genVtxTimeSrc", edm::InputTag("genParticles:t0"));
0597 desc.add<edm::InputTag>("vtxSrc", edm::InputTag("offlinePrimaryVertices4D"));
0598 desc.add<bool>("updateTrackTrajectory", true);
0599 desc.add<bool>("updateTrackExtra", true);
0600 desc.add<bool>("updateTrackHitPattern", true);
0601 desc.add<std::string>("TransientTrackBuilder", "TransientTrackBuilder");
0602 desc.add<std::string>("MTDRecHitBuilder", "MTDRecHitBuilder");
0603 desc.add<std::string>("Propagator", "PropagatorWithMaterialForMTD");
0604 TrackTransformer::fillPSetDescription(transDesc,
0605 false,
0606 "KFFitterForRefitInsideOut",
0607 "KFSmootherForRefitInsideOut",
0608 "PropagatorWithMaterialForMTD",
0609 "alongMomentum",
0610 true,
0611 "WithTrackAngle",
0612 "MuonRecHitBuilder",
0613 "MTDRecHitBuilder");
0614 desc.add<edm::ParameterSetDescription>("TrackTransformer", transDesc);
0615 desc.add<double>("estimatorMaxChi2", 500.);
0616 desc.add<double>("estimatorMaxNSigma", 10.);
0617 desc.add<double>("btlChi2Cut", 50.);
0618 desc.add<double>("btlTimeChi2Cut", 10.);
0619 desc.add<double>("etlChi2Cut", 50.);
0620 desc.add<double>("etlTimeChi2Cut", 10.);
0621 desc.add<bool>("useVertex", false);
0622 desc.add<bool>("useSimVertex", false);
0623 desc.add<double>("dZCut", 0.1);
0624 desc.add<double>("bsTimeSpread", 0.2);
0625 descriptions.add("trackExtenderWithMTDBase", desc);
0626 }
0627
0628 template <class TrackCollection>
0629 template <class H, class T>
0630 void TrackExtenderWithMTDT<TrackCollection>::fillValueMap(edm::Event& iEvent,
0631 const H& handle,
0632 const std::vector<T>& vec,
0633 const edm::EDPutToken& token) const {
0634 auto out = std::make_unique<edm::ValueMap<T>>();
0635 typename edm::ValueMap<T>::Filler filler(*out);
0636 filler.insert(handle, vec.begin(), vec.end());
0637 filler.fill();
0638 iEvent.put(token, std::move(out));
0639 }
0640
0641 template <class TrackCollection>
0642 void TrackExtenderWithMTDT<TrackCollection>::produce(edm::Event& ev, const edm::EventSetup& es) {
0643
0644 Traj2TrackHits t2t;
0645
0646 theTransformer->setServices(es);
0647
0648 TrackingRecHitRefProd hitsRefProd = ev.getRefBeforePut<TrackingRecHitCollection>();
0649 reco::TrackExtraRefProd extrasRefProd = ev.getRefBeforePut<reco::TrackExtraCollection>();
0650
0651 gtg_ = es.getHandle(gtgToken_);
0652
0653 auto geo = es.getTransientHandle(dlgeoToken_);
0654
0655 auto magfield = es.getTransientHandle(magfldToken_);
0656
0657 builder_ = es.getHandle(builderToken_);
0658 hitbuilder_ = es.getHandle(hitbuilderToken_);
0659
0660 auto propH = es.getTransientHandle(propToken_);
0661 const Propagator* prop = propH.product();
0662
0663 auto httopo = es.getTransientHandle(ttopoToken_);
0664 const TrackerTopology& ttopo = *httopo;
0665
0666 auto output = std::make_unique<TrackCollection>();
0667 auto extras = std::make_unique<reco::TrackExtraCollection>();
0668 auto outhits = std::make_unique<edm::OwnVector<TrackingRecHit>>();
0669
0670 std::vector<float> btlMatchChi2;
0671 std::vector<float> etlMatchChi2;
0672 std::vector<float> btlMatchTimeChi2;
0673 std::vector<float> etlMatchTimeChi2;
0674 std::vector<int> npixBarrel;
0675 std::vector<int> npixEndcap;
0676 std::vector<float> pOrigTrkRaw;
0677 std::vector<float> betaOrigTrkRaw;
0678 std::vector<float> t0OrigTrkRaw;
0679 std::vector<float> sigmat0OrigTrkRaw;
0680 std::vector<float> pathLengthsOrigTrkRaw;
0681 std::vector<float> tmtdOrigTrkRaw;
0682 std::vector<float> sigmatmtdOrigTrkRaw;
0683 std::vector<float> tofpiOrigTrkRaw;
0684 std::vector<float> tofkOrigTrkRaw;
0685 std::vector<float> tofpOrigTrkRaw;
0686 std::vector<int> assocOrigTrkRaw;
0687
0688 auto const tracksH = ev.getHandle(tracksToken_);
0689
0690 const auto& trjtrks = ev.get(trajTrackAToken_);
0691
0692
0693 const auto& hits = ev.get(hitsToken_);
0694
0695
0696 const auto& bs = ev.get(bsToken_);
0697
0698 const Vertex* pv = nullptr;
0699 if (useVertex_ && !useSimVertex_) {
0700 auto const& vtxs = ev.get(vtxToken_);
0701 if (!vtxs.empty())
0702 pv = &vtxs[0];
0703 }
0704
0705 std::unique_ptr<math::XYZTLorentzVectorF> genPV(nullptr);
0706 if (useVertex_ && useSimVertex_) {
0707 const auto& genVtxPosition = ev.get(genVtxPositionToken_);
0708 const auto& genVtxTime = ev.get(genVtxTimeToken_);
0709 genPV = std::make_unique<math::XYZTLorentzVectorF>(
0710 genVtxPosition.x(), genVtxPosition.y(), genVtxPosition.z(), genVtxTime);
0711 }
0712
0713 float vtxTime = 0.f;
0714 if (useVertex_) {
0715 if (useSimVertex_ && genPV) {
0716 vtxTime = genPV->t();
0717 } else if (pv)
0718 vtxTime = pv->t();
0719 }
0720
0721 std::vector<unsigned> track_indices;
0722 unsigned itrack = 0;
0723
0724 for (const auto& trjtrk : trjtrks) {
0725 const Trajectory& trajs = *trjtrk.key;
0726 const reco::TrackRef& track = trjtrk.val;
0727
0728 float trackVtxTime = 0.f;
0729 if (useVertex_) {
0730 float dz;
0731 if (useSimVertex_)
0732 dz = std::abs(track->dz(math::XYZPoint(*genPV)));
0733 else
0734 dz = std::abs(track->dz(pv->position()));
0735
0736 if (dz < dzCut_)
0737 trackVtxTime = vtxTime;
0738 }
0739
0740 reco::TransientTrack ttrack(track, magfield.product(), gtg_);
0741 auto thits = theTransformer->getTransientRecHits(ttrack);
0742 TransientTrackingRecHit::ConstRecHitContainer mtdthits;
0743 MTDHitMatchingInfo mBTL, mETL;
0744
0745 if (trajs.isValid()) {
0746
0747 TrajectoryStateOnSurface tsos = builder_->build(track).outermostMeasurementState();
0748 TrajectoryStateClosestToBeamLine tscbl;
0749 bool tscbl_status = getTrajectoryStateClosestToBeamLine(trajs, bs, prop, tscbl);
0750
0751 if (tscbl_status) {
0752 float pmag2 = tscbl.trackStateAtPCA().momentum().mag2();
0753 float pathlength0;
0754 TrackSegments trs0;
0755 trackPathLength(trajs, tscbl, prop, pathlength0, trs0);
0756
0757 const auto& btlhits = tryBTLLayers(tsos,
0758 trajs,
0759 pmag2,
0760 pathlength0,
0761 trs0,
0762 hits,
0763 geo.product(),
0764 magfield.product(),
0765 prop,
0766 bs,
0767 trackVtxTime,
0768 trackVtxTime != 0.f,
0769 mBTL);
0770 mtdthits.insert(mtdthits.end(), btlhits.begin(), btlhits.end());
0771
0772
0773
0774 const auto& etlhits = tryETLLayers(tsos,
0775 trajs,
0776 pmag2,
0777 pathlength0,
0778 trs0,
0779 hits,
0780 geo.product(),
0781 magfield.product(),
0782 prop,
0783 bs,
0784 trackVtxTime,
0785 trackVtxTime != 0.f,
0786 mETL);
0787 mtdthits.insert(mtdthits.end(), etlhits.begin(), etlhits.end());
0788 }
0789 }
0790
0791 auto ordering = checkRecHitsOrdering(thits);
0792 if (ordering == RefitDirection::insideOut) {
0793 thits.insert(thits.end(), mtdthits.begin(), mtdthits.end());
0794 } else {
0795 std::reverse(mtdthits.begin(), mtdthits.end());
0796 mtdthits.insert(mtdthits.end(), thits.begin(), thits.end());
0797 thits.swap(mtdthits);
0798 }
0799
0800 const auto& trajwithmtd =
0801 mtdthits.empty() ? std::vector<Trajectory>(1, trajs) : theTransformer->transform(ttrack, thits);
0802 float pMap = 0.f, betaMap = 0.f, t0Map = 0.f, sigmat0Map = -1.f, pathLengthMap = -1.f, tmtdMap = 0.f,
0803 sigmatmtdMap = -1.f, tofpiMap = 0.f, tofkMap = 0.f, tofpMap = 0.f;
0804 int iMap = -1;
0805
0806 for (const auto& trj : trajwithmtd) {
0807 const auto& thetrj = (updateTraj_ ? trj : trajs);
0808 float pathLength = 0.f, tmtd = 0.f, sigmatmtd = -1.f, tofpi = 0.f, tofk = 0.f, tofp = 0.f;
0809 LogTrace("TrackExtenderWithMTD") << "Refit track " << itrack << " p/pT = " << track->p() << " " << track->pt()
0810 << " eta = " << track->eta();
0811 reco::Track result = buildTrack(track,
0812 thetrj,
0813 trj,
0814 bs,
0815 magfield.product(),
0816 prop,
0817 !trajwithmtd.empty() && !mtdthits.empty(),
0818 pathLength,
0819 tmtd,
0820 sigmatmtd,
0821 tofpi,
0822 tofk,
0823 tofp);
0824 if (result.ndof() >= 0) {
0825
0826 reco::TrackExtra::TrajParams trajParams;
0827 reco::TrackExtra::Chi2sFive chi2s;
0828 size_t hitsstart = outhits->size();
0829 if (updatePattern_) {
0830 t2t(trj, *outhits, trajParams, chi2s);
0831 } else {
0832 t2t(thetrj, *outhits, trajParams, chi2s);
0833 }
0834 size_t hitsend = outhits->size();
0835 extras->push_back(buildTrackExtra(trj));
0836 extras->back().setHits(hitsRefProd, hitsstart, hitsend - hitsstart);
0837 extras->back().setTrajParams(trajParams, chi2s);
0838
0839 output->push_back(result);
0840 btlMatchChi2.push_back(mBTL.hit ? mBTL.estChi2 : -1.f);
0841 etlMatchChi2.push_back(mETL.hit ? mETL.estChi2 : -1.f);
0842 btlMatchTimeChi2.push_back(mBTL.hit ? mBTL.timeChi2 : -1.f);
0843 etlMatchTimeChi2.push_back(mETL.hit ? mETL.timeChi2 : -1.f);
0844 pathLengthMap = pathLength;
0845 tmtdMap = tmtd;
0846 sigmatmtdMap = sigmatmtd;
0847 auto& backtrack = output->back();
0848 iMap = output->size() - 1;
0849 pMap = backtrack.p();
0850 betaMap = backtrack.beta();
0851 t0Map = backtrack.t0();
0852 sigmat0Map = std::copysign(std::sqrt(std::abs(backtrack.covt0t0())), backtrack.covt0t0());
0853 tofpiMap = tofpi;
0854 tofkMap = tofk;
0855 tofpMap = tofp;
0856 reco::TrackExtraRef extraRef(extrasRefProd, extras->size() - 1);
0857 backtrack.setExtra((updateExtra_ ? extraRef : track->extra()));
0858 for (unsigned ihit = hitsstart; ihit < hitsend; ++ihit) {
0859 backtrack.appendHitPattern((*outhits)[ihit], ttopo);
0860 }
0861 npixBarrel.push_back(backtrack.hitPattern().numberOfValidPixelBarrelHits());
0862 npixEndcap.push_back(backtrack.hitPattern().numberOfValidPixelEndcapHits());
0863 LogTrace("TrackExtenderWithMTD") << "tmtd " << tmtdMap << " +/- " << sigmatmtdMap << " t0 " << t0Map << " +/- "
0864 << sigmat0Map << " tof pi/K/p " << tofpiMap << " " << tofkMap << " "
0865 << tofpMap;
0866 } else {
0867 LogTrace("TrackExtenderWithMTD") << "Error in the MTD track refitting. This should not happen";
0868 }
0869 }
0870
0871 pOrigTrkRaw.push_back(pMap);
0872 betaOrigTrkRaw.push_back(betaMap);
0873 t0OrigTrkRaw.push_back(t0Map);
0874 sigmat0OrigTrkRaw.push_back(sigmat0Map);
0875 pathLengthsOrigTrkRaw.push_back(pathLengthMap);
0876 tmtdOrigTrkRaw.push_back(tmtdMap);
0877 sigmatmtdOrigTrkRaw.push_back(sigmatmtdMap);
0878 tofpiOrigTrkRaw.push_back(tofpiMap);
0879 tofkOrigTrkRaw.push_back(tofkMap);
0880 tofpOrigTrkRaw.push_back(tofpMap);
0881 assocOrigTrkRaw.push_back(iMap);
0882
0883 if (iMap == -1) {
0884 btlMatchChi2.push_back(-1.f);
0885 etlMatchChi2.push_back(-1.f);
0886 btlMatchTimeChi2.push_back(-1.f);
0887 etlMatchTimeChi2.push_back(-1.f);
0888 npixBarrel.push_back(-1.f);
0889 npixEndcap.push_back(-1.f);
0890 }
0891
0892 ++itrack;
0893 }
0894
0895 auto outTrksHandle = ev.put(std::move(output));
0896 ev.put(std::move(extras));
0897 ev.put(std::move(outhits));
0898
0899 fillValueMap(ev, tracksH, btlMatchChi2, btlMatchChi2Token_);
0900 fillValueMap(ev, tracksH, etlMatchChi2, etlMatchChi2Token_);
0901 fillValueMap(ev, tracksH, btlMatchTimeChi2, btlMatchTimeChi2Token_);
0902 fillValueMap(ev, tracksH, etlMatchTimeChi2, etlMatchTimeChi2Token_);
0903 fillValueMap(ev, tracksH, npixBarrel, npixBarrelToken_);
0904 fillValueMap(ev, tracksH, npixEndcap, npixEndcapToken_);
0905 fillValueMap(ev, tracksH, pOrigTrkRaw, pOrigTrkToken_);
0906 fillValueMap(ev, tracksH, betaOrigTrkRaw, betaOrigTrkToken_);
0907 fillValueMap(ev, tracksH, t0OrigTrkRaw, t0OrigTrkToken_);
0908 fillValueMap(ev, tracksH, sigmat0OrigTrkRaw, sigmat0OrigTrkToken_);
0909 fillValueMap(ev, tracksH, pathLengthsOrigTrkRaw, pathLengthOrigTrkToken_);
0910 fillValueMap(ev, tracksH, tmtdOrigTrkRaw, tmtdOrigTrkToken_);
0911 fillValueMap(ev, tracksH, sigmatmtdOrigTrkRaw, sigmatmtdOrigTrkToken_);
0912 fillValueMap(ev, tracksH, tofpiOrigTrkRaw, tofpiOrigTrkToken_);
0913 fillValueMap(ev, tracksH, tofkOrigTrkRaw, tofkOrigTrkToken_);
0914 fillValueMap(ev, tracksH, tofpOrigTrkRaw, tofpOrigTrkToken_);
0915 fillValueMap(ev, tracksH, assocOrigTrkRaw, assocOrigTrkToken_);
0916 }
0917
0918 namespace {
0919 bool cmp_for_detset(const unsigned one, const unsigned two) { return one < two; };
0920
0921 void find_hits_in_dets(const MTDTrackingDetSetVector& hits,
0922 const Trajectory& traj,
0923 const DetLayer* layer,
0924 const TrajectoryStateOnSurface& tsos,
0925 const float pmag2,
0926 const float pathlength0,
0927 const TrackSegments& trs0,
0928 const float vtxTime,
0929 const reco::BeamSpot& bs,
0930 const float bsTimeSpread,
0931 const Propagator* prop,
0932 const MeasurementEstimator* estimator,
0933 bool useVtxConstraint,
0934 std::set<MTDHitMatchingInfo>& out) {
0935 pair<bool, TrajectoryStateOnSurface> comp = layer->compatible(tsos, *prop, *estimator);
0936 if (comp.first) {
0937 const vector<DetLayer::DetWithState> compDets = layer->compatibleDets(tsos, *prop, *estimator);
0938 if (!compDets.empty()) {
0939 for (const auto& detWithState : compDets) {
0940 auto range = hits.equal_range(detWithState.first->geographicalId(), cmp_for_detset);
0941 if (range.first == range.second)
0942 continue;
0943
0944 auto pl = prop->propagateWithPath(tsos, detWithState.second.surface());
0945 if (pl.second == 0.)
0946 continue;
0947
0948 const float t_vtx = useVtxConstraint ? vtxTime : 0.f;
0949
0950 constexpr float vtx_res = 0.008f;
0951 const float t_vtx_err = useVtxConstraint ? vtx_res : bsTimeSpread;
0952
0953 float lastpmag2 = trs0.getSegmentPathAndMom2(0).second;
0954
0955 for (auto detitr = range.first; detitr != range.second; ++detitr) {
0956 for (const auto& hit : *detitr) {
0957 auto est = estimator->estimate(detWithState.second, hit);
0958 if (!est.first)
0959 continue;
0960
0961 TrackTofPidInfo tof = computeTrackTofPidInfo(lastpmag2,
0962 std::abs(pl.second),
0963 trs0,
0964 hit.time(),
0965 hit.timeError(),
0966 t_vtx,
0967 t_vtx_err,
0968 false,
0969 TofCalc::kMixd);
0970 MTDHitMatchingInfo mi;
0971 mi.hit = &hit;
0972 mi.estChi2 = est.second;
0973 mi.timeChi2 = tof.dtchi2_best;
0974
0975 out.insert(mi);
0976 }
0977 }
0978 }
0979 }
0980 }
0981 }
0982 }
0983
0984 template <class TrackCollection>
0985 TransientTrackingRecHit::ConstRecHitContainer TrackExtenderWithMTDT<TrackCollection>::tryBTLLayers(
0986 const TrajectoryStateOnSurface& tsos,
0987 const Trajectory& traj,
0988 const float pmag2,
0989 const float pathlength0,
0990 const TrackSegments& trs0,
0991 const MTDTrackingDetSetVector& hits,
0992 const MTDDetLayerGeometry* geo,
0993 const MagneticField* field,
0994 const Propagator* prop,
0995 const reco::BeamSpot& bs,
0996 const float vtxTime,
0997 const bool matchVertex,
0998 MTDHitMatchingInfo& bestHit) const {
0999 const vector<const DetLayer*>& layers = geo->allBTLLayers();
1000
1001 TransientTrackingRecHit::ConstRecHitContainer output;
1002 bestHit = MTDHitMatchingInfo();
1003 for (const DetLayer* ilay : layers)
1004 fillMatchingHits(ilay, tsos, traj, pmag2, pathlength0, trs0, hits, prop, bs, vtxTime, matchVertex, output, bestHit);
1005 return output;
1006 }
1007
1008 template <class TrackCollection>
1009 TransientTrackingRecHit::ConstRecHitContainer TrackExtenderWithMTDT<TrackCollection>::tryETLLayers(
1010 const TrajectoryStateOnSurface& tsos,
1011 const Trajectory& traj,
1012 const float pmag2,
1013 const float pathlength0,
1014 const TrackSegments& trs0,
1015 const MTDTrackingDetSetVector& hits,
1016 const MTDDetLayerGeometry* geo,
1017 const MagneticField* field,
1018 const Propagator* prop,
1019 const reco::BeamSpot& bs,
1020 const float vtxTime,
1021 const bool matchVertex,
1022 MTDHitMatchingInfo& bestHit) const {
1023 const vector<const DetLayer*>& layers = geo->allETLLayers();
1024
1025 TransientTrackingRecHit::ConstRecHitContainer output;
1026 bestHit = MTDHitMatchingInfo();
1027 for (const DetLayer* ilay : layers) {
1028 const BoundDisk& disk = static_cast<const ForwardDetLayer*>(ilay)->specificSurface();
1029 const float diskZ = disk.position().z();
1030
1031 if (tsos.globalPosition().z() * diskZ < 0)
1032 continue;
1033
1034 fillMatchingHits(ilay, tsos, traj, pmag2, pathlength0, trs0, hits, prop, bs, vtxTime, matchVertex, output, bestHit);
1035 }
1036
1037
1038
1039 if (output.size() == 2) {
1040 if (std::abs(output[0]->globalPosition().z()) > std::abs(output[1]->globalPosition().z())) {
1041 std::reverse(output.begin(), output.end());
1042 }
1043 }
1044 return output;
1045 }
1046
1047 template <class TrackCollection>
1048 void TrackExtenderWithMTDT<TrackCollection>::fillMatchingHits(const DetLayer* ilay,
1049 const TrajectoryStateOnSurface& tsos,
1050 const Trajectory& traj,
1051 const float pmag2,
1052 const float pathlength0,
1053 const TrackSegments& trs0,
1054 const MTDTrackingDetSetVector& hits,
1055 const Propagator* prop,
1056 const reco::BeamSpot& bs,
1057 const float& vtxTime,
1058 const bool matchVertex,
1059 TransientTrackingRecHit::ConstRecHitContainer& output,
1060 MTDHitMatchingInfo& bestHit) const {
1061 std::set<MTDHitMatchingInfo> hitsInLayer;
1062 bool hitMatched = false;
1063
1064 using namespace std::placeholders;
1065 auto find_hits = std::bind(find_hits_in_dets,
1066 std::cref(hits),
1067 std::cref(traj),
1068 ilay,
1069 std::cref(tsos),
1070 pmag2,
1071 pathlength0,
1072 trs0,
1073 _1,
1074 std::cref(bs),
1075 bsTimeSpread_,
1076 prop,
1077 theEstimator.get(),
1078 _2,
1079 std::ref(hitsInLayer));
1080
1081 if (useVertex_ && matchVertex)
1082 find_hits(vtxTime, true);
1083 else
1084 find_hits(0, false);
1085
1086 float spaceChi2Cut = ilay->isBarrel() ? btlChi2Cut_ : etlChi2Cut_;
1087 float timeChi2Cut = ilay->isBarrel() ? btlTimeChi2Cut_ : etlTimeChi2Cut_;
1088
1089
1090 if (!hitsInLayer.empty()) {
1091
1092 auto const& firstHit = *hitsInLayer.begin();
1093 if (firstHit.estChi2 < spaceChi2Cut && firstHit.timeChi2 < timeChi2Cut) {
1094 hitMatched = true;
1095 output.push_back(hitbuilder_->build(firstHit.hit));
1096 if (firstHit < bestHit)
1097 bestHit = firstHit;
1098 }
1099 }
1100
1101 if (useVertex_ && matchVertex && !hitMatched) {
1102
1103 hitsInLayer.clear();
1104 find_hits(0, false);
1105 if (!hitsInLayer.empty()) {
1106 auto const& firstHit = *hitsInLayer.begin();
1107 if (firstHit.timeChi2 < timeChi2Cut) {
1108 if (firstHit.estChi2 < spaceChi2Cut) {
1109 hitMatched = true;
1110 output.push_back(hitbuilder_->build(firstHit.hit));
1111 if (firstHit < bestHit)
1112 bestHit = firstHit;
1113 }
1114 }
1115 }
1116 }
1117 }
1118
1119
1120
1121 template <class TrackCollection>
1122 reco::Track TrackExtenderWithMTDT<TrackCollection>::buildTrack(const reco::TrackRef& orig,
1123 const Trajectory& traj,
1124 const Trajectory& trajWithMtd,
1125 const reco::BeamSpot& bs,
1126 const MagneticField* field,
1127 const Propagator* thePropagator,
1128 bool hasMTD,
1129 float& pathLengthOut,
1130 float& tmtdOut,
1131 float& sigmatmtdOut,
1132 float& tofpi,
1133 float& tofk,
1134 float& tofp) const {
1135 TrajectoryStateClosestToBeamLine tscbl;
1136 bool tsbcl_status = getTrajectoryStateClosestToBeamLine(traj, bs, thePropagator, tscbl);
1137
1138 if (!tsbcl_status)
1139 return reco::Track();
1140
1141 GlobalPoint v = tscbl.trackStateAtPCA().position();
1142 math::XYZPoint pos(v.x(), v.y(), v.z());
1143 GlobalVector p = tscbl.trackStateAtPCA().momentum();
1144 math::XYZVector mom(p.x(), p.y(), p.z());
1145
1146 int ndof = traj.ndof();
1147
1148 float t0 = 0.f;
1149 float covt0t0 = -1.f;
1150 pathLengthOut = -1.f;
1151 tmtdOut = 0.f;
1152 sigmatmtdOut = -1.f;
1153 float betaOut = 0.f;
1154 float covbetabeta = -1.f;
1155
1156 auto routput = [&]() {
1157 return reco::Track(traj.chiSquared(),
1158 int(ndof),
1159 pos,
1160 mom,
1161 tscbl.trackStateAtPCA().charge(),
1162 tscbl.trackStateAtPCA().curvilinearError(),
1163 orig->algo(),
1164 reco::TrackBase::undefQuality,
1165 t0,
1166 betaOut,
1167 covt0t0,
1168 covbetabeta);
1169 };
1170
1171
1172 if (hasMTD) {
1173 float pathlength;
1174 TrackSegments trs;
1175 bool validpropagation = trackPathLength(trajWithMtd, bs, thePropagator, pathlength, trs);
1176 float thit = 0.f;
1177 float thiterror = -1.f;
1178 bool validmtd = false;
1179
1180 if (!validpropagation) {
1181 return routput();
1182 }
1183
1184 uint32_t ihitcount(0), ietlcount(0);
1185 for (auto const& hit : trajWithMtd.measurements()) {
1186 if (hit.recHit()->geographicalId().det() == DetId::Forward &&
1187 ForwardSubdetector(hit.recHit()->geographicalId().subdetId()) == FastTime) {
1188 ihitcount++;
1189 if (MTDDetId(hit.recHit()->geographicalId()).mtdSubDetector() == MTDDetId::MTDType::ETL) {
1190 ietlcount++;
1191 }
1192 }
1193 }
1194
1195 auto ihit1 = trajWithMtd.measurements().cbegin();
1196 if (ihitcount == 1) {
1197 const MTDTrackingRecHit* mtdhit = static_cast<const MTDTrackingRecHit*>((*ihit1).recHit()->hit());
1198 thit = mtdhit->time();
1199 thiterror = mtdhit->timeError();
1200 validmtd = true;
1201 } else if (ihitcount == 2 && ietlcount == 2) {
1202 std::pair<float, float> lastStep = trs.getSegmentPathAndMom2(0);
1203 float etlpathlength = std::abs(lastStep.first * c_cm_ns);
1204
1205
1206
1207 if (etlpathlength == 0.f) {
1208 validpropagation = false;
1209 } else {
1210 pathlength -= etlpathlength;
1211 trs.removeFirstSegment();
1212 const MTDTrackingRecHit* mtdhit1 = static_cast<const MTDTrackingRecHit*>((*ihit1).recHit()->hit());
1213 const MTDTrackingRecHit* mtdhit2 = static_cast<const MTDTrackingRecHit*>((*(ihit1 + 1)).recHit()->hit());
1214 TrackTofPidInfo tofInfo = computeTrackTofPidInfo(
1215 lastStep.second, etlpathlength, trs, mtdhit1->time(), mtdhit1->timeError(), 0.f, 0.f, true, TofCalc::kCost);
1216
1217
1218
1219 float err1 = tofInfo.dterror * tofInfo.dterror;
1220 float err2 = mtdhit2->timeError() * mtdhit2->timeError();
1221 if (cms_rounding::roundIfNear0(err1) == 0.f || cms_rounding::roundIfNear0(err2) == 0.f) {
1222 edm::LogError("TrackExtenderWithMTD")
1223 << "MTD tracking hits with zero time uncertainty: " << err1 << " " << err2;
1224 } else {
1225 if ((tofInfo.dt - mtdhit2->time()) * (tofInfo.dt - mtdhit2->time()) < (err1 + err2) * etlTimeChi2Cut_) {
1226
1227
1228
1229
1230 err1 = 1.f / err1;
1231 err2 = 1.f / err2;
1232 thiterror = 1.f / (err1 + err2);
1233 thit = (tofInfo.dt * err1 + mtdhit2->time() * err2) * thiterror;
1234 thiterror = std::sqrt(thiterror);
1235 LogDebug("TrackExtenderWithMTD") << "p trk = " << p.mag() << " ETL hits times/errors: " << mtdhit1->time()
1236 << " +/- " << mtdhit1->timeError() << " , " << mtdhit2->time() << " +/- "
1237 << mtdhit2->timeError() << " extrapolated time1: " << tofInfo.dt << " +/- "
1238 << tofInfo.dterror << " average = " << thit << " +/- " << thiterror;
1239 validmtd = true;
1240 }
1241 }
1242 }
1243 } else {
1244 edm::LogInfo("TrackExtenderWithMTD")
1245 << "MTD hits #" << ihitcount << "ETL hits #" << ietlcount << " anomalous pattern, skipping...";
1246 }
1247
1248 if (validmtd && validpropagation) {
1249
1250 TrackTofPidInfo tofInfo =
1251 computeTrackTofPidInfo(p.mag2(), pathlength, trs, thit, thiterror, 0.f, 0.f, true, TofCalc::kSegm);
1252 pathLengthOut = pathlength;
1253 tmtdOut = thit;
1254 sigmatmtdOut = thiterror;
1255 t0 = tofInfo.dt;
1256 covt0t0 = tofInfo.dterror * tofInfo.dterror;
1257 betaOut = tofInfo.beta_pi;
1258 covbetabeta = tofInfo.betaerror * tofInfo.betaerror;
1259 tofpi = tofInfo.dt_pi;
1260 tofk = tofInfo.dt_k;
1261 tofp = tofInfo.dt_p;
1262 }
1263 }
1264
1265 return routput();
1266 }
1267
1268 template <class TrackCollection>
1269 reco::TrackExtra TrackExtenderWithMTDT<TrackCollection>::buildTrackExtra(const Trajectory& trajectory) const {
1270 static const string metname = "TrackExtenderWithMTD";
1271
1272 const Trajectory::RecHitContainer transRecHits = trajectory.recHits();
1273
1274
1275
1276
1277
1278 TrajectoryStateOnSurface outerTSOS;
1279 TrajectoryStateOnSurface innerTSOS;
1280 unsigned int innerId = 0, outerId = 0;
1281 TrajectoryMeasurement::ConstRecHitPointer outerRecHit;
1282 DetId outerDetId;
1283
1284 if (trajectory.direction() == alongMomentum) {
1285 LogTrace(metname) << "alongMomentum";
1286 outerTSOS = trajectory.lastMeasurement().updatedState();
1287 innerTSOS = trajectory.firstMeasurement().updatedState();
1288 outerId = trajectory.lastMeasurement().recHit()->geographicalId().rawId();
1289 innerId = trajectory.firstMeasurement().recHit()->geographicalId().rawId();
1290 outerRecHit = trajectory.lastMeasurement().recHit();
1291 outerDetId = trajectory.lastMeasurement().recHit()->geographicalId();
1292 } else if (trajectory.direction() == oppositeToMomentum) {
1293 LogTrace(metname) << "oppositeToMomentum";
1294 outerTSOS = trajectory.firstMeasurement().updatedState();
1295 innerTSOS = trajectory.lastMeasurement().updatedState();
1296 outerId = trajectory.firstMeasurement().recHit()->geographicalId().rawId();
1297 innerId = trajectory.lastMeasurement().recHit()->geographicalId().rawId();
1298 outerRecHit = trajectory.firstMeasurement().recHit();
1299 outerDetId = trajectory.firstMeasurement().recHit()->geographicalId();
1300 } else
1301 LogError(metname) << "Wrong propagation direction!";
1302
1303 const GeomDet* outerDet = gtg_->idToDet(outerDetId);
1304 GlobalPoint outerTSOSPos = outerTSOS.globalParameters().position();
1305 bool inside = outerDet->surface().bounds().inside(outerDet->toLocal(outerTSOSPos));
1306
1307 GlobalPoint hitPos =
1308 (outerRecHit->isValid()) ? outerRecHit->globalPosition() : outerTSOS.globalParameters().position();
1309
1310 if (!inside) {
1311 LogTrace(metname) << "The Global Muon outerMostMeasurementState is not compatible with the recHit detector!"
1312 << " Setting outerMost postition to recHit position if recHit isValid: "
1313 << outerRecHit->isValid();
1314 LogTrace(metname) << "From " << outerTSOSPos << " to " << hitPos;
1315 }
1316
1317
1318 GlobalPoint v = (inside) ? outerTSOSPos : hitPos;
1319 GlobalVector p = outerTSOS.globalParameters().momentum();
1320 math::XYZPoint outpos(v.x(), v.y(), v.z());
1321 math::XYZVector outmom(p.x(), p.y(), p.z());
1322
1323 v = innerTSOS.globalParameters().position();
1324 p = innerTSOS.globalParameters().momentum();
1325 math::XYZPoint inpos(v.x(), v.y(), v.z());
1326 math::XYZVector inmom(p.x(), p.y(), p.z());
1327
1328 reco::TrackExtra trackExtra(outpos,
1329 outmom,
1330 true,
1331 inpos,
1332 inmom,
1333 true,
1334 outerTSOS.curvilinearError(),
1335 outerId,
1336 innerTSOS.curvilinearError(),
1337 innerId,
1338 trajectory.direction(),
1339 trajectory.seedRef());
1340
1341 return trackExtra;
1342 }
1343
1344 template <class TrackCollection>
1345 string TrackExtenderWithMTDT<TrackCollection>::dumpLayer(const DetLayer* layer) const {
1346 stringstream output;
1347
1348 const BoundSurface* sur = nullptr;
1349 const BoundCylinder* bc = nullptr;
1350 const BoundDisk* bd = nullptr;
1351
1352 sur = &(layer->surface());
1353 if ((bc = dynamic_cast<const BoundCylinder*>(sur))) {
1354 output << " Cylinder of radius: " << bc->radius() << endl;
1355 } else if ((bd = dynamic_cast<const BoundDisk*>(sur))) {
1356 output << " Disk at: " << bd->position().z() << endl;
1357 }
1358 return output.str();
1359 }
1360
1361
1362 #include <FWCore/Framework/interface/MakerMacros.h>
1363 #include "DataFormats/TrackReco/interface/Track.h"
1364 #include "DataFormats/TrackReco/interface/TrackFwd.h"
1365 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
1366 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
1367 typedef TrackExtenderWithMTDT<reco::TrackCollection> TrackExtenderWithMTD;
1368
1369 DEFINE_FWK_MODULE(TrackExtenderWithMTD);