Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-02-21 23:14:04 

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