TrackAssociatorProducers/plugins/TrackAssociatorByPositionImpl.cc

0001 #include "TrackAssociatorByPositionImpl.h"
0002 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0003
0004 #include <TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h>
0005 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
0006
0007 #include "DataFormats/Math/interface/deltaR.h"
0008 //#include "PhysicsTools/Utilities/interface/DeltaR.h"
0009
0010 using namespace edm;
0011 using namespace reco;
0012
0013 TrajectoryStateOnSurface TrackAssociatorByPositionImpl::getState(const TrackingParticleRef& st,
0014                                                                  const SimHitTPAssociationList& simHitsTPAssoc) const {
0015   using SimHitTPPair = std::pair<TrackingParticleRef, TrackPSimHitRef>;
0016   SimHitTPPair clusterTPpairWithDummyTP(st, TrackPSimHitRef());  //SimHit is dummy:
0017   // sorting only the cluster is needed
0018   auto range = std::equal_range(
0019       simHitsTPAssoc.begin(),
0020       simHitsTPAssoc.end(),
0021       clusterTPpairWithDummyTP,
0022       [](const SimHitTPPair& iLHS, const SimHitTPPair& iRHS) -> bool { return iLHS.first.key() > iRHS.first.key(); });
0023
0024   //  TrackingParticle* simtrack = const_cast<TrackingParticle*>(&st);
0025   //loop over PSimHits
0026   const PSimHit* psimhit = nullptr;
0027   const BoundPlane* plane = nullptr;
0028   double dLim = thePositionMinimumDistance;
0029
0030   //    look for the further most hit beyond a certain limit
0031   auto start = range.first;
0032   auto end = range.second;
0033   LogDebug("TrackAssociatorByPositionImpl") << range.second - range.first << " PSimHits.";
0034
0035   unsigned int count = 0;
0036   for (auto ip = start; ip != end; ++ip) {
0037     TrackPSimHitRef psit = ip->second;
0038
0039     //get the detid
0040     DetId dd(psit->detUnitId());
0041
0042     if (!theConsiderAllSimHits && dd.det() != DetId::Tracker)
0043       continue;
0044
0045     LogDebug("TrackAssociatorByPositionImpl") << count++ << "] PSimHit on: " << dd.rawId();
0046     //get the surface from the global geometry
0047     const GeomDet* gd = theGeometry->idToDet(dd);
0048     if (!gd) {
0049       edm::LogError("TrackAssociatorByPositionImpl") << "no geomdet for: " << dd.rawId() << ". will skip.";
0050       continue;
0051     }
0052     double d = gd->surface().toGlobal(psit->localPosition()).mag();
0053     if (d > dLim) {
0054       dLim = d;
0055       psimhit = &(*psit);
0056       plane = &gd->surface();
0057     }
0058   }
0059
0060   if (psimhit && plane) {
0061     //build a trajectorystate on this surface
0062     SurfaceSideDefinition::SurfaceSide surfaceside = SurfaceSideDefinition::atCenterOfSurface;
0063     GlobalPoint initialPoint = plane->toGlobal(psimhit->localPosition());
0064     GlobalVector initialMomentum = plane->toGlobal(psimhit->momentumAtEntry());
0065     int initialCharge = (psimhit->particleType() > 0) ? -1 : 1;
0066     CartesianTrajectoryError initialCartesianErrors;  //no error at initial state
0067     const GlobalTrajectoryParameters initialParameters(
0068         initialPoint, initialMomentum, initialCharge, thePropagator->magneticField());
0069     return TrajectoryStateOnSurface(initialParameters, initialCartesianErrors, *plane, surfaceside);
0070   } else {
0071     //    edm::LogError("TrackAssociatorByPositionImpl")<<"no corresponding PSimHit for a tracking particle. will fail.";
0072     return TrajectoryStateOnSurface();
0073   }
0074 }
0075
0076 FreeTrajectoryState TrackAssociatorByPositionImpl::getState(const reco::Track& track) const {
0077   //may be you want to do more than that if track does not go to IP
0078   return trajectoryStateTransform::initialFreeState(track, thePropagator->magneticField());
0079 }
0080
0081 double TrackAssociatorByPositionImpl::quality(const TrajectoryStateOnSurface& tr,
0082                                               const TrajectoryStateOnSurface& sim) const {
0083   switch (theMethod) {
0084     case Method::chi2: {
0085       AlgebraicVector5 v(tr.localParameters().vector() - sim.localParameters().vector());
0086       AlgebraicSymMatrix55 m(tr.localError().matrix());
0087       int ierr = !m.Invert();
0088       if (ierr != 0)
0089         edm::LogInfo("TrackAssociatorByPositionImpl") << "error inverting the error matrix:\n" << m;
0090       double est = ROOT::Math::Similarity(v, m);
0091       return est;
0092       break;
0093     }
0094     case Method::dist: {
0095       return (tr.globalPosition() - sim.globalPosition()).mag();
0096       break;
0097     }
0098     case Method::momdr: {
0099       return (deltaR<double>(tr.globalDirection().eta(),
0100                              tr.globalDirection().phi(),
0101                              sim.globalDirection().eta(),
0102                              sim.globalDirection().phi()));
0103       break;
0104     }
0105     case Method::posdr: {
0106       return (deltaR<double>(
0107           tr.globalPosition().eta(), tr.globalPosition().phi(), sim.globalPosition().eta(), sim.globalPosition().phi()));
0108       break;
0109     }
0110   }
0111   //should never be reached
0112   edm::LogError("TrackAssociatorByPositionImpl")
0113       << "option: " << static_cast<int>(theMethod) << " has not been recognized. association has no meaning.";
0114   return -1;
0115 }
0116
0117 RecoToSimCollection TrackAssociatorByPositionImpl::associateRecoToSim(
0118     const edm::RefToBaseVector<reco::Track>& tCH, const edm::RefVector<TrackingParticleCollection>& tPCH) const {
0119   RecoToSimCollection outputCollection(productGetter_);
0120   //for each reco track find a matching tracking particle
0121   std::pair<unsigned int, unsigned int> minPair;
0122   const double dQmin_default = 1542543;
0123   double dQmin = dQmin_default;
0124
0125   //cdj edm::Handle<SimHitTPAssociationProducer::SimHitTPAssociationList> simHitsTPAssoc;
0126
0127   //warning: make sure the TP collection used in the map is the same used in the associator!
0128   //e->getByLabel(_simHitTpMapTag,simHitsTPAssoc);
0129
0130   for (unsigned int Ti = 0; Ti != tCH.size(); ++Ti) {
0131     //initial state (initial OR inner OR outter)
0132     FreeTrajectoryState iState = getState(*(tCH)[Ti]);
0133
0134     bool atLeastOne = false;
0135     //    for each tracking particle, find a state position and the plane to propagate the track to.
0136     for (unsigned int TPi = 0; TPi != tPCH.size(); ++TPi) {
0137       //get a state in the muon system
0138       TrajectoryStateOnSurface simReferenceState = getState((tPCH)[TPi], *theSimHitsTPAssoc);
0139       if (!simReferenceState.isValid())
0140         continue;
0141
0142       //propagate the TRACK to the surface
0143       TrajectoryStateOnSurface trackReferenceState = thePropagator->propagate(iState, simReferenceState.surface());
0144       if (!trackReferenceState.isValid())
0145         continue;
0146
0147       //comparison
0148       double dQ = quality(trackReferenceState, simReferenceState);
0149       if (dQ < theQCut) {
0150         atLeastOne = true;
0151         outputCollection.insert(tCH[Ti],
0152                                 std::make_pair(tPCH[TPi], -dQ));  //association map with quality, is order greater-first
0153         edm::LogVerbatim("TrackAssociatorByPositionImpl")
0154             << "track number: " << Ti << " associated with dQ: " << dQ << " to TrackingParticle number: " << TPi;
0155       }
0156       if (dQ < dQmin) {
0157         dQmin = dQ;
0158         minPair = std::make_pair(Ti, TPi);
0159       }
0160     }  //loop over tracking particles
0161     if (theMinIfNoMatch && !atLeastOne && dQmin != dQmin_default) {
0162       outputCollection.insert(tCH[minPair.first], std::make_pair(tPCH[minPair.second], -dQmin));
0163     }
0164   }  //loop over tracks
0165   outputCollection.post_insert();
0166   return outputCollection;
0167 }
0168
0169 SimToRecoCollection TrackAssociatorByPositionImpl::associateSimToReco(
0170     const edm::RefToBaseVector<reco::Track>& tCH, const edm::RefVector<TrackingParticleCollection>& tPCH) const {
0171   SimToRecoCollection outputCollection(productGetter_);
0172   //for each tracking particle, find matching tracks.
0173
0174   std::pair<unsigned int, unsigned int> minPair;
0175   const double dQmin_default = 1542543;
0176   double dQmin = dQmin_default;
0177
0178   //edm::Handle<SimHitTPAssociationProducer::SimHitTPAssociationList> simHitsTPAssoc;
0179
0180   //warning: make sure the TP collection used in the map is the same used in the associator!
0181   //e->getByLabel(_simHitTpMapTag,simHitsTPAssoc);
0182
0183   for (unsigned int TPi = 0; TPi != tPCH.size(); ++TPi) {
0184     //get a state in the muon system
0185     TrajectoryStateOnSurface simReferenceState = getState((tPCH)[TPi], *theSimHitsTPAssoc);
0186
0187     if (!simReferenceState.isValid())
0188       continue;
0189     bool atLeastOne = false;
0190     //  propagate every track from any state (initial, inner, outter) to the surface
0191     //  and make the position test
0192     for (unsigned int Ti = 0; Ti != tCH.size(); ++Ti) {
0193       //initial state
0194       FreeTrajectoryState iState = getState(*(tCH)[Ti]);
0195
0196       //propagation to surface
0197       TrajectoryStateOnSurface trackReferenceState = thePropagator->propagate(iState, simReferenceState.surface());
0198       if (!trackReferenceState.isValid())
0199         continue;
0200
0201       //comparison
0202       double dQ = quality(trackReferenceState, simReferenceState);
0203       if (dQ < theQCut) {
0204         atLeastOne = true;
0205         outputCollection.insert(tPCH[TPi],
0206                                 std::make_pair(tCH[Ti], -dQ));  //association map with quality, is order greater-first
0207         edm::LogVerbatim("TrackAssociatorByPositionImpl")
0208             << "TrackingParticle number: " << TPi << " associated with dQ: " << dQ << " to track number: " << Ti;
0209       }
0210       if (dQ < dQmin) {
0211         dQmin = dQ;
0212         minPair = std::make_pair(TPi, Ti);
0213       }
0214     }  //loop over tracks
0215     if (theMinIfNoMatch && !atLeastOne && dQmin != dQmin_default) {
0216       outputCollection.insert(tPCH[minPair.first], std::make_pair(tCH[minPair.second], -dQmin));
0217     }
0218   }  //loop over tracking particles
0219
0220   outputCollection.post_insert();
0221   return outputCollection;
0222 }