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