Back to home page

Project CMSSW displayed by LXR
 
 

    

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);  // [mm/ns] -> [cm/ns]
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     //Operator used to sort the hits while performing the matching step at the MTD
0080     inline bool operator<(const MTDHitMatchingInfo& m2) const {
0081       //only for good matching in time use estChi2, otherwise use mostly time compatibility
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);  // observed upper limit on nSegments
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       // remove previously calculated sigmaTofs
0146       sigmaTofs_.clear();
0147 
0148       // compute sigma(tof) on each segment first by propagating sigma(p)
0149       // also add diagonal terms to sigmatof
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       // compute sigma on sum of tofs assuming full correlation between segments
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           // sigma(t) = sigma(p) * |dt/dp| = sigma(p) * DeltaL/c * m^2 / (p^2 * E)
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;  //assume by default the pi hypothesis
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       // only add sigma(TOF) if not considering mass hp. uncertainty
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       //*TODO* deal with heavier nucleons and/or BSM case here?
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     // get the state closest to the beamline
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     //add pathlength layer by layer
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       // sigma(p) from curvilinear error (on q/p)
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     //add distance from bs to first measurement
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 }  // namespace
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       // maybe perp2?
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   //this produces pieces of the track extra
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   //MTD hits DetSet
0804   const auto& hits = ev.get(hitsToken_);
0805 
0806   //beam spot
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();  //already in ns
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       // get the outermost trajectory point on the track
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         // in the future this should include an intermediate refit before propagating to the ETL
0890         // for now it is ok
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         /// setup the track extras
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);  // this fills the output hit collection
0961         } else {
0962           t2t(thetrj, *outhits, trajParams, chi2s);
0963         }
0964         size_t hitsend = outhits->size();
0965         extras->push_back(buildTrackExtra(trj));  // always push back the fully built extra, update by setting in track
0966         extras->back().setHits(hitsRefProd, hitsstart, hitsend - hitsstart);
0967         extras->back().setTrajParams(trajParams, chi2s);
0968         //create the track
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());  // save R of the outermost hit for BTL, z for ETL.
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,  //put vtx error by hand for the moment
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;  //use the chi2 for the best matching hypothesis
1160 
1161               out.insert(mi);
1162             }
1163           }
1164         }
1165       }
1166     }
1167   }
1168 }  // namespace
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;  // only propagate to the disk that's on the same side
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   // the ETL hits order must be from the innermost to the outermost
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   //just take the first hit because the hits are sorted on their matching quality
1283   if (!hitsInLayer.empty()) {
1284     //check hits to pass minimum quality matching requirements
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     //try a second search with beamspot hypothesis
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 //below is unfortunately ripped from other places but
1326 //since track producer doesn't know about MTD we have to do this
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;  // if there is no MTD flag the pathlength with -1
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   //compute path length for time backpropagation, using first MTD hit for the momentum
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       // The information of the two ETL hits is combined and attributed to the innermost hit
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         // Protect against incompatible times
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             // Subtract the ETL time of flight from the outermost measurement, and combine it in a weighted average with the innermost
1443             // the mass ambiguity related uncertainty on the time of flight is added as an additional uncertainty
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             // if back extrapolated time of the outermost measurement not compatible with the innermost, keep the one with smallest error
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       //here add the PID uncertainty for later use in the 1st step of 4D vtx reconstruction
1479       TrackTofPidInfo tofInfo = computeTrackTofPidInfo(
1480           p.mag2(), pathlength, trs, thit, thiterror, 0.f, 0.f, true, TofCalc::kSegm, SigmaTofCalc::kCost);
1481 
1482       pathLengthOut = pathlength;  // set path length if we've got a timing hit
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   // put the collection of TrackingRecHit in the event
1509 
1510   // sets the outermost and innermost TSOSs
1511   // ToDo: validation for track states with MTD
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   //build the TrackExtra
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 //define this as a plug-in
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);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!