File indexing completed on 2024-04-06 12:26:52
0001
0002
0003
0004
0005
0006
0007 #include "RecoMuon/CosmicMuonProducer/interface/GlobalCosmicMuonTrajectoryBuilder.h"
0008 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011
0012 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
0013 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0014
0015 #include "DataFormats/TrackReco/interface/Track.h"
0016 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0017 #include "DataFormats/TrackReco/interface/TrackExtraFwd.h"
0018 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0019 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
0020 #include "TrackingTools/GeomPropagators/interface/StateOnTrackerBound.h"
0021 #include "DataFormats/Math/interface/deltaPhi.h"
0022
0023 using namespace std;
0024 using namespace edm;
0025
0026
0027
0028
0029 GlobalCosmicMuonTrajectoryBuilder::GlobalCosmicMuonTrajectoryBuilder(const edm::ParameterSet& par,
0030 const MuonServiceProxy* service,
0031 edm::ConsumesCollector& iC)
0032 : theService(service),
0033 theTrackerRecHitBuilderToken(
0034 iC.esConsumes(edm::ESInputTag("", par.getParameter<string>("TrackerRecHitBuilder")))),
0035 theMuonRecHitBuilderToken(iC.esConsumes(edm::ESInputTag("", par.getParameter<string>("MuonRecHitBuilder")))) {
0036 ParameterSet smootherPSet = par.getParameter<ParameterSet>("SmootherParameters");
0037 theSmoother = new CosmicMuonSmoother(smootherPSet, theService);
0038
0039 ParameterSet trackMatcherPSet = par.getParameter<ParameterSet>("GlobalMuonTrackMatcher");
0040 theTrackMatcher = new GlobalMuonTrackMatcher(trackMatcherPSet, theService);
0041
0042 theTkTrackToken = iC.consumes<reco::TrackCollection>(par.getParameter<InputTag>("TkTrackCollectionLabel"));
0043
0044 thePropagatorName = par.getParameter<string>("Propagator");
0045 category_ = "Muon|RecoMuon|CosmicMuon|GlobalCosmicMuonTrajectoryBuilder";
0046 }
0047
0048
0049
0050
0051
0052 GlobalCosmicMuonTrajectoryBuilder::~GlobalCosmicMuonTrajectoryBuilder() {
0053 if (theSmoother)
0054 delete theSmoother;
0055 if (theTrackMatcher)
0056 delete theTrackMatcher;
0057 }
0058
0059
0060
0061
0062 void GlobalCosmicMuonTrajectoryBuilder::setEvent(const edm::Event& event) {
0063 event.getByToken(theTkTrackToken, theTrackerTracks);
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077 theTrackerRecHitBuilder = theService->eventSetup().getHandle(theTrackerRecHitBuilderToken);
0078 theMuonRecHitBuilder = theService->eventSetup().getHandle(theMuonRecHitBuilderToken);
0079 }
0080
0081
0082
0083
0084 MuonCandidate::CandidateContainer GlobalCosmicMuonTrajectoryBuilder::trajectories(const TrackCand& muCand) {
0085 MuonCandidate::CandidateContainer result;
0086
0087 if (!theTrackerTracks.isValid()) {
0088 LogTrace(category_) << "Tracker Track collection is invalid!!!";
0089 return result;
0090 }
0091
0092 LogTrace(category_) << "Found " << theTrackerTracks->size() << " tracker Tracks";
0093 if (theTrackerTracks->empty())
0094 return result;
0095
0096 vector<TrackCand> matched = match(muCand, theTrackerTracks);
0097
0098 LogTrace(category_) << "TrackMatcher found " << matched.size() << "tracker tracks matched";
0099
0100 if (matched.empty())
0101 return result;
0102 reco::TrackRef tkTrack = matched.front().second;
0103
0104 if (tkTrack.isNull())
0105 return result;
0106 reco::TrackRef muTrack = muCand.second;
0107
0108 ConstRecHitContainer muRecHits;
0109
0110 if (muCand.first == nullptr || !muCand.first->isValid()) {
0111 muRecHits = getTransientRecHits(*muTrack);
0112 } else {
0113 muRecHits = muCand.first->recHits();
0114 }
0115
0116 LogTrace(category_) << "mu RecHits: " << muRecHits.size();
0117
0118 ConstRecHitContainer tkRecHits = getTransientRecHits(*tkTrack);
0119
0120
0121
0122
0123
0124
0125
0126 ConstRecHitContainer hits;
0127 LogTrace(category_) << "tk RecHits: " << tkRecHits.size();
0128
0129
0130
0131
0132 sortHits(hits, muRecHits, tkRecHits);
0133
0134
0135
0136
0137
0138 TrajectoryStateOnSurface muonState1 = trajectoryStateTransform::innerStateOnSurface(
0139 *muTrack, *theService->trackingGeometry(), &*theService->magneticField());
0140 TrajectoryStateOnSurface tkState1 = trajectoryStateTransform::innerStateOnSurface(
0141 *tkTrack, *theService->trackingGeometry(), &*theService->magneticField());
0142
0143 TrajectoryStateOnSurface muonState2 = trajectoryStateTransform::outerStateOnSurface(
0144 *muTrack, *theService->trackingGeometry(), &*theService->magneticField());
0145 TrajectoryStateOnSurface tkState2 = trajectoryStateTransform::outerStateOnSurface(
0146 *tkTrack, *theService->trackingGeometry(), &*theService->magneticField());
0147
0148 TrajectoryStateOnSurface firstState1 =
0149 (muonState1.globalPosition().y() > tkState1.globalPosition().y()) ? muonState1 : tkState1;
0150 TrajectoryStateOnSurface firstState2 =
0151 (muonState2.globalPosition().y() > tkState2.globalPosition().y()) ? muonState2 : tkState2;
0152
0153 TrajectoryStateOnSurface firstState =
0154 (firstState1.globalPosition().y() > firstState2.globalPosition().y()) ? firstState1 : firstState2;
0155
0156 GlobalPoint front, back;
0157 if (!hits.empty()) {
0158 front = hits.front()->globalPosition();
0159 back = hits.back()->globalPosition();
0160 if ((front.perp() < 130 && fabs(front.z()) < 300) || (back.perp() < 130 && fabs(back.z()) < 300)) {
0161 if (hits.front()->globalPosition().perp() > hits.back()->globalPosition().perp())
0162 reverse(hits.begin(), hits.end());
0163 tkState1 = trajectoryStateTransform::innerStateOnSurface(
0164 *tkTrack, *theService->trackingGeometry(), &*theService->magneticField());
0165 tkState2 = trajectoryStateTransform::outerStateOnSurface(
0166 *tkTrack, *theService->trackingGeometry(), &*theService->magneticField());
0167 firstState = (tkState1.globalPosition().perp() < tkState2.globalPosition().perp()) ? tkState1 : tkState2;
0168 }
0169 }
0170 if (!firstState.isValid())
0171 return result;
0172 LogTrace(category_) << "firstTSOS pos: " << firstState.globalPosition() << "mom: " << firstState.globalMomentum();
0173
0174
0175
0176 TrajectorySeed seed;
0177 vector<Trajectory> refitted = theSmoother->trajectories(seed, hits, firstState);
0178
0179 if (refitted.empty()) {
0180 LogTrace(category_) << "smoothing trajectories fail";
0181 refitted = theSmoother->fit(seed, hits, firstState);
0182 }
0183
0184 if (refitted.empty()) {
0185 LogTrace(category_) << "refit fail";
0186 return result;
0187 }
0188
0189 auto myTraj = std::make_unique<Trajectory>(refitted.front());
0190
0191 const std::vector<TrajectoryMeasurement>& mytms = myTraj->measurements();
0192 LogTrace(category_) << "measurements in final trajectory " << mytms.size();
0193 LogTrace(category_) << "Orignally there are " << tkTrack->found() << " tk rhs and " << muTrack->found() << " mu rhs.";
0194
0195 if (mytms.size() <= tkTrack->found()) {
0196 LogTrace(category_) << "insufficient measurements. skip... ";
0197 return result;
0198 }
0199
0200 result.push_back(std::make_unique<MuonCandidate>(std::move(myTraj), muTrack, tkTrack));
0201 LogTrace(category_) << "final global cosmic muon: ";
0202 for (std::vector<TrajectoryMeasurement>::const_iterator itm = mytms.begin(); itm != mytms.end(); ++itm) {
0203 LogTrace(category_) << "updated pos " << itm->updatedState().globalPosition() << "mom "
0204 << itm->updatedState().globalMomentum();
0205 }
0206 return result;
0207 }
0208
0209 void GlobalCosmicMuonTrajectoryBuilder::sortHits(ConstRecHitContainer& hits,
0210 ConstRecHitContainer& muonHits,
0211 ConstRecHitContainer& tkHits) {
0212 if (tkHits.empty()) {
0213 LogTrace(category_) << "No valid tracker hits";
0214 return;
0215 }
0216 if (muonHits.empty()) {
0217 LogTrace(category_) << "No valid muon hits";
0218 return;
0219 }
0220
0221 ConstRecHitContainer::const_iterator frontTkHit = tkHits.begin();
0222 ConstRecHitContainer::const_iterator backTkHit = tkHits.end() - 1;
0223 while (!(*frontTkHit)->isValid() && frontTkHit != backTkHit) {
0224 frontTkHit++;
0225 }
0226 while (!(*backTkHit)->isValid() && backTkHit != frontTkHit) {
0227 backTkHit--;
0228 }
0229
0230 ConstRecHitContainer::const_iterator frontMuHit = muonHits.begin();
0231 ConstRecHitContainer::const_iterator backMuHit = muonHits.end() - 1;
0232 while (!(*frontMuHit)->isValid() && frontMuHit != backMuHit) {
0233 frontMuHit++;
0234 }
0235 while (!(*backMuHit)->isValid() && backMuHit != frontMuHit) {
0236 backMuHit--;
0237 }
0238
0239 if (frontTkHit == backTkHit) {
0240 LogTrace(category_) << "No valid tracker hits";
0241 return;
0242 }
0243 if (frontMuHit == backMuHit) {
0244 LogTrace(category_) << "No valid muon hits";
0245 return;
0246 }
0247
0248 GlobalPoint frontTkPos = (*frontTkHit)->globalPosition();
0249 GlobalPoint backTkPos = (*backTkHit)->globalPosition();
0250
0251 GlobalPoint frontMuPos = (*frontMuHit)->globalPosition();
0252 GlobalPoint backMuPos = (*backMuHit)->globalPosition();
0253
0254
0255 if (frontTkPos.y() < backTkPos.y()) {
0256 reverse(tkHits.begin(), tkHits.end());
0257 }
0258
0259 if (frontMuPos.y() < backMuPos.y()) {
0260 reverse(muonHits.begin(), muonHits.end());
0261 }
0262
0263
0264
0265
0266
0267
0268
0269
0270
0271
0272 ConstRecHitContainer::iterator middlepoint = muonHits.begin();
0273 bool insertInMiddle = false;
0274
0275 for (ConstRecHitContainer::iterator ihit = muonHits.begin(); ihit != muonHits.end() - 1; ihit++) {
0276 GlobalPoint ipos = (*ihit)->globalPosition();
0277 GlobalPoint nextpos = (*(ihit + 1))->globalPosition();
0278 if ((ipos - nextpos).mag() < 100.0)
0279 continue;
0280
0281 GlobalPoint middle((ipos.x() + nextpos.x()) / 2, (ipos.y() + nextpos.y()) / 2, (ipos.z() + nextpos.z()) / 2);
0282 LogTrace(category_) << "ipos " << ipos << "nextpos" << nextpos << " middle " << middle << endl;
0283 if ((middle.perp() < ipos.perp()) && (middle.perp() < nextpos.perp())) {
0284 LogTrace(category_) << "found middlepoint" << endl;
0285 middlepoint = ihit;
0286 insertInMiddle = true;
0287 break;
0288 }
0289 }
0290
0291
0292 if (insertInMiddle) {
0293 GlobalPoint jointpointpos = (*middlepoint)->globalPosition();
0294 LogTrace(category_) << "jointpoint " << jointpointpos << endl;
0295 if ((frontTkPos - jointpointpos).mag() >
0296 (backTkPos - jointpointpos).mag()) {
0297 reverse(tkHits.begin(), tkHits.end());
0298 }
0299 muonHits.insert(middlepoint + 1, tkHits.begin(), tkHits.end());
0300 hits = muonHits;
0301 } else {
0302 if (frontTkPos.y() < frontMuPos.y()) {
0303 LogTrace(category_) << "insert at the end " << frontTkPos << frontMuPos << endl;
0304
0305 hits = muonHits;
0306 hits.insert(hits.end(), tkHits.begin(), tkHits.end());
0307 } else {
0308 LogTrace(category_) << "insert at the beginning " << frontTkPos << frontMuPos << endl;
0309 hits = tkHits;
0310 hits.insert(hits.end(), muonHits.begin(), muonHits.end());
0311 }
0312 }
0313 }
0314
0315 #include "RecoTracker/TransientTrackingRecHit/interface/TkTransientTrackingRecHitBuilder.h"
0316
0317 TransientTrackingRecHit::ConstRecHitContainer GlobalCosmicMuonTrajectoryBuilder::getTransientRecHits(
0318 const reco::Track& track) const {
0319 TransientTrackingRecHit::ConstRecHitContainer result;
0320
0321 auto hitCloner = static_cast<TkTransientTrackingRecHitBuilder const*>(theTrackerRecHitBuilder.product())->cloner();
0322 TrajectoryStateOnSurface currTsos = trajectoryStateTransform::innerStateOnSurface(
0323 track, *theService->trackingGeometry(), &*theService->magneticField());
0324 for (trackingRecHit_iterator hit = track.recHitsBegin(); hit != track.recHitsEnd(); ++hit) {
0325 if ((*hit)->isValid()) {
0326 DetId recoid = (*hit)->geographicalId();
0327 if (recoid.det() == DetId::Tracker) {
0328 TrajectoryStateOnSurface predTsos =
0329 theService->propagator(thePropagatorName)
0330 ->propagate(currTsos, theService->trackingGeometry()->idToDet(recoid)->surface());
0331 LogTrace(category_) << "predtsos " << predTsos.isValid();
0332 if (predTsos.isValid()) {
0333 currTsos = predTsos;
0334 result.emplace_back(hitCloner(**hit, predTsos));
0335 }
0336 } else if (recoid.det() == DetId::Muon) {
0337 result.push_back(theMuonRecHitBuilder->build(&**hit));
0338 }
0339 }
0340 }
0341 return result;
0342 }
0343
0344 std::vector<GlobalCosmicMuonTrajectoryBuilder::TrackCand> GlobalCosmicMuonTrajectoryBuilder::match(
0345 const TrackCand& mu, const edm::Handle<reco::TrackCollection>& tktracks) {
0346 std::vector<TrackCand> result;
0347
0348 TrajectoryStateOnSurface innerTsos = trajectoryStateTransform::innerStateOnSurface(
0349 *(mu.second), *theService->trackingGeometry(), &*theService->magneticField());
0350 TrajectoryStateOnSurface outerTsos = trajectoryStateTransform::outerStateOnSurface(
0351 *(mu.second), *theService->trackingGeometry(), &*theService->magneticField());
0352
0353 vector<TrackCand> tkTrackCands;
0354 for (reco::TrackCollection::size_type i = 0; i < theTrackerTracks->size(); ++i) {
0355 reco::TrackRef tkTrack(theTrackerTracks, i);
0356 TrackCand tkCand = TrackCand((Trajectory*)nullptr, tkTrack);
0357 tkTrackCands.push_back(tkCand);
0358 LogTrace(category_) << "chisq is " << theTrackMatcher->match(mu, tkCand, 0, 0);
0359 LogTrace(category_) << "d is " << theTrackMatcher->match(mu, tkCand, 1, 0);
0360 LogTrace(category_) << "r_pos is " << theTrackMatcher->match(mu, tkCand, 2, 0);
0361 }
0362
0363
0364 if (tkTrackCands.size() <= 1) {
0365 return tkTrackCands;
0366 }
0367
0368
0369
0370
0371 GlobalPoint innerPos = innerTsos.globalPosition();
0372 GlobalPoint outerPos = outerTsos.globalPosition();
0373
0374 if ((innerPos.basicVector().dot(innerTsos.globalMomentum().basicVector()) *
0375 outerPos.basicVector().dot(outerTsos.globalMomentum().basicVector()) >
0376 0)) {
0377 GlobalPoint geoInnerPos = (innerPos.mag() < outerPos.mag()) ? innerPos : outerPos;
0378 LogTrace(category_) << "geoInnerPos Mu " << geoInnerPos << endl;
0379
0380
0381
0382
0383 for (vector<TrackCand>::const_iterator itkCand = tkTrackCands.begin(); itkCand != tkTrackCands.end(); ++itkCand) {
0384 reco::TrackRef tkTrack = itkCand->second;
0385
0386 GlobalPoint tkInnerPos(tkTrack->innerPosition().x(), tkTrack->innerPosition().y(), tkTrack->innerPosition().z());
0387 GlobalPoint tkOuterPos(tkTrack->outerPosition().x(), tkTrack->outerPosition().y(), tkTrack->outerPosition().z());
0388 LogTrace(category_) << "tkTrack " << tkInnerPos << " " << tkOuterPos << endl;
0389
0390 float closetDistance11 = (geoInnerPos - tkInnerPos).mag();
0391 float closetDistance12 = (geoInnerPos - tkOuterPos).mag();
0392 float closetDistance1 = (closetDistance11 < closetDistance12) ? closetDistance11 : closetDistance12;
0393 LogTrace(category_) << "closetDistance1 " << closetDistance1 << endl;
0394
0395 if (true || !isTraversing(*tkTrack)) {
0396 bool keep = true;
0397 for (vector<TrackCand>::const_iterator itkCand2 = tkTrackCands.begin(); itkCand2 != tkTrackCands.end();
0398 ++itkCand2) {
0399 if (itkCand2 == itkCand)
0400 continue;
0401 reco::TrackRef tkTrack2 = itkCand2->second;
0402
0403 GlobalPoint tkInnerPos2(
0404 tkTrack2->innerPosition().x(), tkTrack2->innerPosition().y(), tkTrack2->innerPosition().z());
0405 GlobalPoint tkOuterPos2(
0406 tkTrack2->outerPosition().x(), tkTrack2->outerPosition().y(), tkTrack2->outerPosition().z());
0407 LogTrace(category_) << "tkTrack2 " << tkInnerPos2 << " " << tkOuterPos2 << endl;
0408
0409 float farthestDistance21 = (geoInnerPos - tkInnerPos2).mag();
0410 float farthestDistance22 = (geoInnerPos - tkOuterPos2).mag();
0411 float farthestDistance2 = (farthestDistance21 > farthestDistance22) ? farthestDistance21 : farthestDistance22;
0412 LogTrace(category_) << "farthestDistance2 " << farthestDistance2 << endl;
0413
0414 if (closetDistance1 > farthestDistance2 - 1e-3) {
0415 keep = false;
0416 break;
0417 }
0418 }
0419 if (keep)
0420 result.push_back(*itkCand);
0421 else
0422 LogTrace(category_) << "The Track is on different hemisphere" << endl;
0423 } else {
0424 result.push_back(*itkCand);
0425 }
0426 }
0427 if (result.empty()) {
0428
0429 result = tkTrackCands;
0430 }
0431 } else {
0432 result = tkTrackCands;
0433 }
0434
0435
0436 vector<TrackCand> matched_trackerTracks = theTrackMatcher->match(mu, result);
0437
0438 LogTrace(category_) << "TrackMatcher found " << matched_trackerTracks.size() << "tracker tracks matched";
0439
0440
0441 if (matched_trackerTracks.size() < 2) {
0442 return matched_trackerTracks;
0443 } else {
0444
0445
0446
0447 result.clear();
0448 TrackCand bestMatch;
0449
0450 double quality = 1e6;
0451 double max_quality = 1e6;
0452 for (vector<TrackCand>::const_iterator iter = matched_trackerTracks.begin(); iter != matched_trackerTracks.end();
0453 iter++) {
0454 quality = theTrackMatcher->match(mu, *iter, 1, 0);
0455 LogTrace(category_) << " quality of tracker track is " << quality;
0456 if (quality < max_quality) {
0457 max_quality = quality;
0458 bestMatch = (*iter);
0459 }
0460 }
0461 LogTrace(category_) << " Picked tracker track with quality " << max_quality;
0462 result.push_back(bestMatch);
0463 return result;
0464 }
0465 }
0466
0467 bool GlobalCosmicMuonTrajectoryBuilder::isTraversing(const reco::Track& track) const {
0468 trackingRecHit_iterator firstValid;
0469 for (trackingRecHit_iterator hit = track.recHitsBegin(); hit != track.recHitsEnd(); ++hit) {
0470 if ((*hit)->isValid()) {
0471 firstValid = hit;
0472 break;
0473 }
0474 }
0475
0476 trackingRecHit_iterator lastValid;
0477 for (trackingRecHit_iterator hit = track.recHitsEnd() - 1; hit != track.recHitsBegin() - 1; --hit) {
0478 if ((*hit)->isValid()) {
0479 lastValid = hit;
0480 break;
0481 }
0482 }
0483
0484 GlobalPoint posFirst = theService->trackingGeometry()->idToDet((*firstValid)->geographicalId())->position();
0485
0486 GlobalPoint posLast = theService->trackingGeometry()->idToDet((*lastValid)->geographicalId())->position();
0487
0488 GlobalPoint middle(
0489 (posFirst.x() + posLast.x()) / 2, (posFirst.y() + posLast.y()) / 2, (posFirst.z() + posLast.z()) / 2);
0490
0491 if ((middle.mag() < posFirst.mag()) && (middle.mag() < posLast.mag())) {
0492 return true;
0493 }
0494 return false;
0495 }