File indexing completed on 2023-03-17 11:20:25
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include "RecoMuon/GlobalTrackingTools/interface/GlobalTrajectoryBuilderBase.h"
0021
0022
0023
0024
0025
0026 #include <iostream>
0027 #include <algorithm>
0028 #include <limits>
0029
0030
0031
0032
0033
0034 #include "FWCore/Framework/interface/Event.h"
0035 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0036 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0037
0038 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0039 #include "TrackingTools/TrackFitters/interface/RecHitLessByDet.h"
0040 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0041 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0042
0043 #include "DataFormats/Math/interface/deltaR.h"
0044
0045 #include "DataFormats/DetId/interface/DetId.h"
0046 #include "DataFormats/MuonDetId/interface/DTChamberId.h"
0047 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0048 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
0049
0050 #include "DataFormats/TrackReco/interface/Track.h"
0051 #include "DataFormats/TrackReco/interface/TrackExtraFwd.h"
0052 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0053 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0054
0055 #include "RecoMuon/GlobalTrackingTools/interface/GlobalMuonTrackMatcher.h"
0056 #include "RecoMuon/GlobalTrackingTools/interface/GlobalMuonRefitter.h"
0057 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHitBuilder.h"
0058 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"
0059 #include "RecoMuon/TrackingTools/interface/MuonCandidate.h"
0060 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0061
0062 #include "RecoMuon/GlobalTrackingTools/interface/MuonTrackingRegionBuilder.h"
0063 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0064 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0065 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0066
0067 #include "RecoTracker/TkTrackingRegions/interface/TkTrackingRegionsMargin.h"
0068 #include "RecoTracker/TkMSParametrization/interface/PixelRecoRange.h"
0069
0070 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0071 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0072
0073
0074
0075
0076 GlobalTrajectoryBuilderBase::GlobalTrajectoryBuilderBase(const edm::ParameterSet& par,
0077 const MuonServiceProxy* service,
0078 edm::ConsumesCollector& iC)
0079 : theLayerMeasurements(nullptr),
0080 theTrackTransformer(nullptr),
0081 theRegionBuilder(nullptr),
0082 theService(service),
0083 theGlbRefitter(nullptr) {
0084 theCategory =
0085 par.getUntrackedParameter<std::string>("Category", "Muon|RecoMuon|GlobalMuon|GlobalTrajectoryBuilderBase");
0086
0087 edm::ParameterSet trackMatcherPSet = par.getParameter<edm::ParameterSet>("GlobalMuonTrackMatcher");
0088 theTrackMatcher = new GlobalMuonTrackMatcher(trackMatcherPSet, theService);
0089
0090 theTrackerPropagatorName = par.getParameter<std::string>("TrackerPropagator");
0091
0092 edm::ParameterSet trackTransformerPSet = par.getParameter<edm::ParameterSet>("TrackTransformer");
0093 theTrackTransformer = new TrackTransformer(trackTransformerPSet, iC);
0094
0095 edm::ParameterSet regionBuilderPSet = par.getParameter<edm::ParameterSet>("MuonTrackingRegionBuilder");
0096
0097 theRegionBuilder = new MuonTrackingRegionBuilder(regionBuilderPSet, iC);
0098
0099
0100 edm::ParameterSet refitterParameters = par.getParameter<edm::ParameterSet>("GlbRefitterParameters");
0101 theGlbRefitter = new GlobalMuonRefitter(refitterParameters, theService, iC);
0102
0103 theMuonHitsOption = refitterParameters.getParameter<int>("MuonHitsOption");
0104 theRefitFlag = refitterParameters.getParameter<bool>("RefitFlag");
0105
0106 theTrackerRecHitBuilderToken =
0107 iC.esConsumes(edm::ESInputTag("", par.getParameter<std::string>("TrackerRecHitBuilder")));
0108 theMuonRecHitBuilderToken = iC.esConsumes(edm::ESInputTag("", par.getParameter<std::string>("MuonRecHitBuilder")));
0109 theTopoToken = iC.esConsumes();
0110 theFieldToken = iC.esConsumes();
0111 theMSMakerToken = iC.esConsumes();
0112
0113 theRPCInTheFit = par.getParameter<bool>("RefitRPCHits");
0114
0115 theTECxScale = par.getParameter<double>("ScaleTECxFactor");
0116 theTECyScale = par.getParameter<double>("ScaleTECyFactor");
0117 thePtCut = par.getParameter<double>("PtCut");
0118 thePCut = par.getParameter<double>("PCut");
0119 }
0120
0121
0122
0123
0124 GlobalTrajectoryBuilderBase::~GlobalTrajectoryBuilderBase() {
0125 if (theTrackMatcher)
0126 delete theTrackMatcher;
0127 if (theRegionBuilder)
0128 delete theRegionBuilder;
0129 if (theTrackTransformer)
0130 delete theTrackTransformer;
0131 if (theGlbRefitter)
0132 delete theGlbRefitter;
0133 }
0134
0135
0136
0137
0138 void GlobalTrajectoryBuilderBase::setEvent(const edm::Event& event) {
0139 theEvent = &event;
0140
0141 theTrackTransformer->setServices(theService->eventSetup());
0142 theRegionBuilder->setEvent(event, theService->eventSetup());
0143
0144 theGlbRefitter->setEvent(event);
0145 theGlbRefitter->setServices(theService->eventSetup());
0146
0147 theTrackerRecHitBuilder = &theService->eventSetup().getData(theTrackerRecHitBuilderToken);
0148 theMuonRecHitBuilder = &theService->eventSetup().getData(theMuonRecHitBuilderToken);
0149
0150
0151 theTopo = &theService->eventSetup().getData(theTopoToken);
0152
0153 theField = &theService->eventSetup().getData(theFieldToken);
0154 theMSMaker = &theService->eventSetup().getData(theMSMakerToken);
0155 }
0156
0157
0158
0159
0160 MuonCandidate::CandidateContainer GlobalTrajectoryBuilderBase::build(const TrackCand& staCand,
0161 MuonCandidate::CandidateContainer& tkTrajs) const {
0162 LogTrace(theCategory) << " Begin Build" << std::endl;
0163
0164
0165 if (tkTrajs.empty())
0166 return CandidateContainer();
0167
0168
0169 CandidateContainer refittedResult;
0170 ConstRecHitContainer muonRecHits = getTransientRecHits(*(staCand.second));
0171
0172
0173 if ((muonRecHits.size() > 1) &&
0174 (muonRecHits.front()->globalPosition().mag() > muonRecHits.back()->globalPosition().mag())) {
0175 LogTrace(theCategory) << " reverse order: ";
0176 }
0177
0178 for (auto&& it : tkTrajs) {
0179
0180 LogTrace(theCategory) << " Track p and pT " << it->trackerTrack()->p() << " " << it->trackerTrack()->pt();
0181 if (it->trackerTrack()->p() < thePCut || it->trackerTrack()->pt() < thePtCut)
0182 continue;
0183
0184
0185 if (theRefitFlag) {
0186 ConstRecHitContainer trackerRecHits;
0187 if (it->trackerTrack().isNonnull()) {
0188 trackerRecHits = getTransientRecHits(*it->trackerTrack());
0189 } else {
0190 LogDebug(theCategory) << " NEED HITS FROM TRAJ";
0191 }
0192
0193
0194
0195 if (std::abs(it->trackerTrack()->eta()) > 0.95 && std::abs(it->trackerTrack()->eta()) < 1.15 &&
0196 it->trackerTrack()->pt() < 60) {
0197 if (theTECxScale < 0 || theTECyScale < 0)
0198 trackerRecHits = selectTrackerHits(trackerRecHits);
0199 else
0200 fixTEC(trackerRecHits, theTECxScale, theTECyScale);
0201 }
0202
0203 RefitDirection recHitDir = checkRecHitsOrdering(trackerRecHits);
0204 if (recHitDir == outToIn)
0205 reverse(trackerRecHits.begin(), trackerRecHits.end());
0206
0207 reco::TransientTrack tTT(it->trackerTrack(), &*theService->magneticField(), theService->trackingGeometry());
0208 TrajectoryStateOnSurface innerTsos = tTT.innermostMeasurementState();
0209
0210 edm::RefToBase<TrajectorySeed> tmpSeed;
0211 if (it->trackerTrack()->seedRef().isAvailable())
0212 tmpSeed = it->trackerTrack()->seedRef();
0213
0214 if (!innerTsos.isValid()) {
0215 LogTrace(theCategory) << " inner Trajectory State is invalid. ";
0216 continue;
0217 }
0218
0219 innerTsos.rescaleError(100.);
0220
0221 TC refitted0, refitted1;
0222 std::unique_ptr<Trajectory> tkTrajectory;
0223
0224
0225 if (!(it->trackerTrajectory() && it->trackerTrajectory()->isValid())) {
0226 refitted0 = theTrackTransformer->transform(it->trackerTrack());
0227 if (!refitted0.empty())
0228 tkTrajectory = std::make_unique<Trajectory>(*(refitted0.begin()));
0229 else
0230 LogDebug(theCategory) << " Failed to load tracker track trajectory";
0231 } else
0232 tkTrajectory = it->releaseTrackerTrajectory();
0233 if (tkTrajectory)
0234 tkTrajectory->setSeedRef(tmpSeed);
0235
0236
0237 ConstRecHitContainer allRecHits = trackerRecHits;
0238 allRecHits.insert(allRecHits.end(), muonRecHits.begin(), muonRecHits.end());
0239 refitted1 = theGlbRefitter->refit(*it->trackerTrack(), tTT, allRecHits, theMuonHitsOption, theTopo);
0240 LogTrace(theCategory) << " This track-sta refitted to " << refitted1.size() << " trajectories";
0241
0242 std::unique_ptr<Trajectory> glbTrajectory1;
0243 if (!refitted1.empty())
0244 glbTrajectory1 = std::make_unique<Trajectory>(*(refitted1.begin()));
0245 else
0246 LogDebug(theCategory) << " Failed to load global track trajectory 1";
0247 if (glbTrajectory1)
0248 glbTrajectory1->setSeedRef(tmpSeed);
0249
0250 if (glbTrajectory1 && tkTrajectory) {
0251 refittedResult.emplace_back(std::make_unique<MuonCandidate>(
0252 std::move(glbTrajectory1), it->muonTrack(), it->trackerTrack(), std::move(tkTrajectory)));
0253 }
0254 } else {
0255 edm::RefToBase<TrajectorySeed> tmpSeed;
0256 if (it->trackerTrack()->seedRef().isAvailable())
0257 tmpSeed = it->trackerTrack()->seedRef();
0258
0259 TC refitted0;
0260 std::unique_ptr<Trajectory> tkTrajectory;
0261 if (!(it->trackerTrajectory() && it->trackerTrajectory()->isValid())) {
0262 refitted0 = theTrackTransformer->transform(it->trackerTrack());
0263 if (!refitted0.empty()) {
0264 tkTrajectory = std::make_unique<Trajectory>(*(refitted0.begin()));
0265 } else
0266 LogDebug(theCategory) << " Failed to load tracker track trajectory";
0267 } else
0268 tkTrajectory = it->releaseTrackerTrajectory();
0269 std::unique_ptr<Trajectory> cpy;
0270 if (tkTrajectory) {
0271 tkTrajectory->setSeedRef(tmpSeed);
0272 cpy = std::make_unique<Trajectory>(*tkTrajectory);
0273 }
0274
0275 refittedResult.emplace_back(std::make_unique<MuonCandidate>(
0276 std::move(tkTrajectory), it->muonTrack(), it->trackerTrack(), std::move(cpy)));
0277 }
0278 }
0279
0280
0281 CandidateContainer selectedResult;
0282 std::unique_ptr<MuonCandidate> tmpCand;
0283 double minProb = std::numeric_limits<double>::max();
0284
0285 for (auto&& cand : refittedResult) {
0286 double prob = trackProbability(*cand->trajectory());
0287 LogTrace(theCategory) << " refitted-track-sta with pT " << cand->trackerTrack()->pt() << " has probability "
0288 << prob;
0289
0290 if (prob < minProb or not tmpCand) {
0291 minProb = prob;
0292 tmpCand = std::move(cand);
0293 }
0294 }
0295
0296 if (tmpCand)
0297 selectedResult.push_back(std::move(tmpCand));
0298
0299 refittedResult.clear();
0300
0301 return selectedResult;
0302 }
0303
0304
0305
0306
0307 std::vector<GlobalTrajectoryBuilderBase::TrackCand> GlobalTrajectoryBuilderBase::chooseRegionalTrackerTracks(
0308 const TrackCand& staCand, const std::vector<TrackCand>& tkTs) {
0309
0310 RectangularEtaPhiTrackingRegion regionOfInterest = defineRegionOfInterest(staCand.second);
0311
0312
0313
0314
0315
0316 std::vector<TrackCand> result;
0317
0318 double deltaR_max = 1.0;
0319
0320 for (auto&& is : tkTs) {
0321 double deltaR_tmp = deltaR(static_cast<double>(regionOfInterest.direction().eta()),
0322 static_cast<double>(regionOfInterest.direction().phi()),
0323 is.second->eta(),
0324 is.second->phi());
0325
0326
0327
0328 if (deltaR_tmp < deltaR_max) {
0329 TrackCand tmpCand = TrackCand(is);
0330 result.push_back(tmpCand);
0331 }
0332 }
0333
0334 return result;
0335 }
0336
0337
0338
0339
0340 RectangularEtaPhiTrackingRegion GlobalTrajectoryBuilderBase::defineRegionOfInterest(
0341 const reco::TrackRef& staTrack) const {
0342 std::unique_ptr<RectangularEtaPhiTrackingRegion> region1 = theRegionBuilder->region(staTrack);
0343
0344 TkTrackingRegionsMargin<float> etaMargin(std::abs(region1->etaRange().min() - region1->etaRange().mean()),
0345 std::abs(region1->etaRange().max() - region1->etaRange().mean()));
0346
0347 RectangularEtaPhiTrackingRegion region2(region1->direction(),
0348 region1->origin(),
0349 region1->ptMin(),
0350 region1->originRBound(),
0351 region1->originZBound(),
0352 etaMargin,
0353 region1->phiMargin(),
0354 *theField,
0355 theMSMaker);
0356
0357 return region2;
0358 }
0359
0360
0361
0362
0363 double GlobalTrajectoryBuilderBase::trackProbability(const Trajectory& track) const {
0364 if (track.ndof() > 0 && track.chiSquared() > 0) {
0365 return -LnChiSquaredProbability(track.chiSquared(), track.ndof());
0366 } else {
0367 return 0.0;
0368 }
0369 }
0370
0371
0372
0373
0374 void GlobalTrajectoryBuilderBase::printHits(const ConstRecHitContainer& hits) const {
0375 LogTrace(theCategory) << "Used RecHits: " << hits.size();
0376 for (auto&& ir : hits) {
0377 if (!ir->isValid()) {
0378 LogTrace(theCategory) << "invalid RecHit";
0379 continue;
0380 }
0381
0382 const GlobalPoint& pos = ir->globalPosition();
0383
0384 LogTrace(theCategory) << "r = " << sqrt(pos.x() * pos.x() + pos.y() * pos.y()) << " z = " << pos.z()
0385 << " dimension = " << ir->dimension() << " " << ir->det()->geographicalId().det() << " "
0386 << ir->det()->subDetector();
0387 }
0388 }
0389
0390
0391
0392
0393 GlobalTrajectoryBuilderBase::RefitDirection GlobalTrajectoryBuilderBase::checkRecHitsOrdering(
0394 const TransientTrackingRecHit::ConstRecHitContainer& recHits) const {
0395 if (!recHits.empty()) {
0396 ConstRecHitContainer::const_iterator frontHit = recHits.begin();
0397 ConstRecHitContainer::const_iterator backHit = recHits.end() - 1;
0398 while (!(*frontHit)->isValid() && frontHit != backHit) {
0399 frontHit++;
0400 }
0401 while (!(*backHit)->isValid() && backHit != frontHit) {
0402 backHit--;
0403 }
0404
0405 double rFirst = (*frontHit)->globalPosition().mag();
0406 double rLast = (*backHit)->globalPosition().mag();
0407
0408 if (rFirst < rLast)
0409 return inToOut;
0410 else if (rFirst > rLast)
0411 return outToIn;
0412 else {
0413 edm::LogError(theCategory) << "Impossible to determine the rechits order" << std::endl;
0414 return undetermined;
0415 }
0416 } else {
0417 edm::LogError(theCategory) << "Impossible to determine the rechits order" << std::endl;
0418 return undetermined;
0419 }
0420 }
0421
0422
0423
0424
0425 GlobalTrajectoryBuilderBase::ConstRecHitContainer GlobalTrajectoryBuilderBase::selectTrackerHits(
0426 const ConstRecHitContainer& all) const {
0427 int nTEC(0);
0428
0429 ConstRecHitContainer hits;
0430 for (auto&& i : all) {
0431 if (!i->isValid())
0432 continue;
0433 if (i->det()->geographicalId().det() == DetId::Tracker &&
0434 i->det()->geographicalId().subdetId() == StripSubdetector::TEC) {
0435 nTEC++;
0436 } else {
0437 hits.push_back(i);
0438 }
0439 if (nTEC > 1)
0440 return all;
0441 }
0442
0443 return hits;
0444 }
0445
0446
0447
0448
0449 void GlobalTrajectoryBuilderBase::fixTEC(ConstRecHitContainer& all, double scl_x, double scl_y) const {
0450 int nTEC(0);
0451 ConstRecHitContainer::iterator lone_tec;
0452
0453 for (ConstRecHitContainer::iterator i = all.begin(); i != all.end(); i++) {
0454 if (!(*i)->isValid())
0455 continue;
0456
0457 if ((*i)->det()->geographicalId().det() == DetId::Tracker &&
0458 (*i)->det()->geographicalId().subdetId() == StripSubdetector::TEC) {
0459 lone_tec = i;
0460 nTEC++;
0461
0462 if ((i + 1) != all.end() && (*(i + 1))->isValid() &&
0463 (*(i + 1))->det()->geographicalId().det() == DetId::Tracker &&
0464 (*(i + 1))->det()->geographicalId().subdetId() == StripSubdetector::TEC) {
0465 nTEC++;
0466 break;
0467 }
0468 }
0469
0470 if (nTEC > 1)
0471 break;
0472 }
0473
0474 int hitDet = (*lone_tec)->hit()->geographicalId().det();
0475 int hitSubDet = (*lone_tec)->hit()->geographicalId().subdetId();
0476 if (nTEC == 1 && (*lone_tec)->hit()->isValid() && hitDet == DetId::Tracker && hitSubDet == StripSubdetector::TEC) {
0477
0478 const SiStripRecHit2D* strip = dynamic_cast<const SiStripRecHit2D*>((*lone_tec)->hit());
0479 if (strip && strip->det()) {
0480 LocalPoint pos = strip->localPosition();
0481 if ((*lone_tec)->detUnit()) {
0482 const StripTopology* topology = dynamic_cast<const StripTopology*>(&(*lone_tec)->detUnit()->topology());
0483 if (topology) {
0484
0485 float angle = topology->stripAngle(topology->strip((*lone_tec)->hit()->localPosition()));
0486 LocalError error = strip->localPositionError();
0487 LocalError rotError = error.rotate(angle);
0488 LocalError scaledError(rotError.xx() * scl_x * scl_x, 0, rotError.yy() * scl_y * scl_y);
0489 error = scaledError.rotate(-angle);
0490
0491
0492
0493
0494 SiStripRecHit2D* st = new SiStripRecHit2D(pos, error, *strip->det(), strip->cluster());
0495 *lone_tec = MuonTransientTrackingRecHit::build((*lone_tec)->det(), st);
0496 }
0497 }
0498 }
0499 }
0500 }
0501
0502 #include "RecoTracker/TransientTrackingRecHit/interface/TkTransientTrackingRecHitBuilder.h"
0503
0504
0505
0506 TransientTrackingRecHit::ConstRecHitContainer GlobalTrajectoryBuilderBase::getTransientRecHits(
0507 const reco::Track& track) const {
0508 TransientTrackingRecHit::ConstRecHitContainer result;
0509
0510 TrajectoryStateOnSurface currTsos = trajectoryStateTransform::innerStateOnSurface(
0511 track, *theService->trackingGeometry(), &*theService->magneticField());
0512
0513 auto tkbuilder = static_cast<TkTransientTrackingRecHitBuilder const*>(theTrackerRecHitBuilder);
0514 auto hitCloner = tkbuilder->cloner();
0515 for (trackingRecHit_iterator hit = track.recHitsBegin(); hit != track.recHitsEnd(); ++hit) {
0516 if ((*hit)->isValid()) {
0517 DetId recoid = (*hit)->geographicalId();
0518 if (recoid.det() == DetId::Tracker) {
0519 if (!(*hit)->hasPositionAndError()) {
0520 TrajectoryStateOnSurface predTsos =
0521 theService->propagator(theTrackerPropagatorName)
0522 ->propagate(currTsos, theService->trackingGeometry()->idToDet(recoid)->surface());
0523
0524 if (!predTsos.isValid()) {
0525 edm::LogError("MissingTransientHit")
0526 << "Could not get a tsos on the hit surface. We will miss a tracking hit.";
0527 continue;
0528 }
0529 currTsos = predTsos;
0530 auto h = (**hit).cloneForFit(*tkbuilder->geometry()->idToDet((**hit).geographicalId()));
0531 result.emplace_back(hitCloner.makeShared(h, predTsos));
0532 } else {
0533 result.push_back((*hit)->cloneSH());
0534 }
0535 } else if (recoid.det() == DetId::Muon) {
0536 if ((*hit)->geographicalId().subdetId() == 3 && !theRPCInTheFit) {
0537 LogDebug(theCategory) << "RPC Rec Hit discarded";
0538 continue;
0539 }
0540 result.push_back(theMuonRecHitBuilder->build(&**hit));
0541 }
0542 }
0543 }
0544
0545 return result;
0546 }