Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-06-17 01:30:44 

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