File indexing completed on 2024-04-06 12:26:52
0001 #include "RecoMuon/CosmicMuonProducer/interface/CosmicMuonTrajectoryBuilder.h"
0002
0003
0004
0005
0006
0007
0008
0009
0010 #include "FWCore/Framework/interface/Frameworkfwd.h"
0011 #include "FWCore/Framework/interface/MakerMacros.h"
0012
0013
0014 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0015 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0016 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0017 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
0018 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
0019 #include "FWCore/Framework/interface/Event.h"
0020 #include "FWCore/Framework/interface/EventSetup.h"
0021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0022 #include "FWCore/Utilities/interface/InputTag.h"
0023 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0024 #include "RecoMuon/Navigation/interface/DirectMuonNavigation.h"
0025 #include "RecoMuon/MeasurementDet/interface/MuonDetLayerMeasurements.h"
0026 #include "RecoMuon/TrackingTools/interface/MuonBestMeasurementFinder.h"
0027 #include "RecoMuon/TrackingTools/interface/MuonTrajectoryUpdator.h"
0028 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0029 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0030 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"
0031 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0032 #include "DataFormats/GeometrySurface/interface/PlaneBuilder.h"
0033 #include "DataFormats/GeometryCommonDetAlgo/interface/PerpendicularBoundPlaneBuilder.h"
0034 #include "RecoMuon/Records/interface/MuonRecoGeometryRecord.h"
0035 #include "DataFormats/CSCRecHit/interface/CSCRecHit2D.h"
0036 #include "DataFormats/DTRecHit/interface/DTRecHit1D.h"
0037
0038 #include <algorithm>
0039
0040 using namespace edm;
0041 using namespace std;
0042
0043 CosmicMuonTrajectoryBuilder::CosmicMuonTrajectoryBuilder(const edm::ParameterSet& par,
0044 const MuonServiceProxy* service,
0045 edm::ConsumesCollector& iC)
0046 : theService(service) {
0047 thePropagatorName = par.getParameter<string>("Propagator");
0048
0049 bool enableDTMeasurement = par.getParameter<bool>("EnableDTMeasurement");
0050 bool enableCSCMeasurement = par.getParameter<bool>("EnableCSCMeasurement");
0051 bool enableRPCMeasurement = par.getParameter<bool>("EnableRPCMeasurement");
0052
0053
0054 InputTag DTRecSegmentLabel = par.getParameter<InputTag>("DTRecSegmentLabel");
0055
0056
0057 InputTag CSCRecSegmentLabel = par.getParameter<InputTag>("CSCRecSegmentLabel");
0058
0059
0060 InputTag RPCRecSegmentLabel = par.getParameter<InputTag>("RPCRecSegmentLabel");
0061
0062
0063
0064
0065 theLayerMeasurements = new MuonDetLayerMeasurements(DTRecSegmentLabel,
0066 CSCRecSegmentLabel,
0067 RPCRecSegmentLabel,
0068 edm::InputTag(),
0069 edm::InputTag(),
0070 iC,
0071 enableDTMeasurement,
0072 enableCSCMeasurement,
0073 enableRPCMeasurement,
0074 false,
0075 false);
0076
0077 ParameterSet muonUpdatorPSet = par.getParameter<ParameterSet>("MuonTrajectoryUpdatorParameters");
0078
0079 theNavigation = nullptr;
0080 theUpdator = new MuonTrajectoryUpdator(muonUpdatorPSet, insideOut);
0081
0082 theBestMeasurementFinder = new MuonBestMeasurementFinder();
0083
0084 ParameterSet muonBackwardUpdatorPSet = par.getParameter<ParameterSet>("BackwardMuonTrajectoryUpdatorParameters");
0085
0086 theBKUpdator = new MuonTrajectoryUpdator(muonBackwardUpdatorPSet, outsideIn);
0087
0088 theTraversingMuonFlag = par.getParameter<bool>("BuildTraversingMuon");
0089
0090 theStrict1LegFlag = par.getParameter<bool>("Strict1Leg");
0091
0092 ParameterSet smootherPSet = par.getParameter<ParameterSet>("MuonSmootherParameters");
0093
0094 theNavigationPSet = par.getParameter<ParameterSet>("MuonNavigationParameters");
0095
0096 theSmoother = new CosmicMuonSmoother(smootherPSet, theService);
0097
0098 theNTraversing = 0;
0099 theNSuccess = 0;
0100 theCacheId_DG = 0;
0101 category_ = "Muon|RecoMuon|CosmicMuon|CosmicMuonTrajectoryBuilder";
0102 }
0103
0104 CosmicMuonTrajectoryBuilder::~CosmicMuonTrajectoryBuilder() {
0105 LogTrace(category_) << "CosmicMuonTrajectoryBuilder dtor called";
0106 if (theUpdator)
0107 delete theUpdator;
0108 if (theBKUpdator)
0109 delete theBKUpdator;
0110 if (theLayerMeasurements)
0111 delete theLayerMeasurements;
0112 if (theSmoother)
0113 delete theSmoother;
0114 if (theNavigation)
0115 delete theNavigation;
0116 delete theBestMeasurementFinder;
0117
0118 LogTrace(category_) << "CosmicMuonTrajectoryBuilder Traversing: " << theNSuccess << "/" << theNTraversing;
0119 }
0120
0121 void CosmicMuonTrajectoryBuilder::setEvent(const edm::Event& event) {
0122 theLayerMeasurements->setEvent(event);
0123
0124
0125 unsigned long long newCacheId_DG = theService->eventSetup().get<MuonRecoGeometryRecord>().cacheIdentifier();
0126 if (newCacheId_DG != theCacheId_DG) {
0127 LogTrace(category_) << "Muon Reco Geometry changed!";
0128 theCacheId_DG = newCacheId_DG;
0129 if (theNavigation)
0130 delete theNavigation;
0131 theNavigation = new DirectMuonNavigation(theService->detLayerGeometry(), theNavigationPSet);
0132 }
0133 }
0134
0135 MuonTrajectoryBuilder::TrajectoryContainer CosmicMuonTrajectoryBuilder::trajectories(const TrajectorySeed& seed) {
0136 TrajectoryContainer emptyContainer;
0137
0138 MuonPatternRecoDumper debug;
0139
0140 PTrajectoryStateOnDet ptsd1(seed.startingState());
0141 DetId did(ptsd1.detId());
0142 const BoundPlane& bp = theService->trackingGeometry()->idToDet(did)->surface();
0143 TrajectoryStateOnSurface lastTsos =
0144 trajectoryStateTransform::transientState(ptsd1, &bp, &*theService->magneticField());
0145 LogTrace(category_) << "Seed: mom " << lastTsos.globalMomentum() << "pos: " << lastTsos.globalPosition();
0146 LogTrace(category_) << "Seed: mom eta " << lastTsos.globalMomentum().eta()
0147 << "pos eta: " << lastTsos.globalPosition().eta();
0148
0149 bool beamhaloFlag = ((did.subdetId() == MuonSubdetId::CSC) && fabs(lastTsos.globalMomentum().eta()) > 4.0);
0150
0151 vector<const DetLayer*> navLayers;
0152
0153 if (did.subdetId() == MuonSubdetId::DT) {
0154
0155 navLayers = navigation()->compatibleLayers(*(lastTsos.freeState()), alongMomentum);
0156 } else if (beamhaloFlag || (theTraversingMuonFlag && theStrict1LegFlag)) {
0157
0158 navLayers = navigation()->compatibleEndcapLayers(*(lastTsos.freeState()), alongMomentum);
0159 } else {
0160 navLayers = navigation()->compatibleLayers(*(lastTsos.freeState()), alongMomentum);
0161 }
0162
0163 LogTrace(category_) << "found " << navLayers.size() << " compatible DetLayers for the Seed";
0164
0165 if (navLayers.empty())
0166 return emptyContainer;
0167
0168 vector<DetWithState> detsWithStates;
0169 LogTrace(category_) << "Compatible layers: ";
0170 for (vector<const DetLayer*>::const_iterator layer = navLayers.begin(); layer != navLayers.end(); layer++) {
0171 LogTrace(category_) << debug.dumpMuonId((*layer)->basicComponents().front()->geographicalId())
0172 << debug.dumpLayer(*layer);
0173 }
0174
0175 detsWithStates = navLayers.front()->compatibleDets(lastTsos, *propagator(), *(updator()->estimator()));
0176 LogTrace(category_) << "Number of compatible dets: " << detsWithStates.size() << endl;
0177
0178 if (!detsWithStates.empty()) {
0179
0180 if (detsWithStates.front().second.isValid()) {
0181 LogTrace(category_) << "New starting TSOS is on det: " << endl;
0182 LogTrace(category_) << debug.dumpMuonId(detsWithStates.front().first->geographicalId())
0183 << debug.dumpLayer(navLayers.front());
0184 lastTsos = detsWithStates.front().second;
0185 LogTrace(category_) << "Seed after extrapolation: mom " << lastTsos.globalMomentum()
0186 << "pos: " << lastTsos.globalPosition();
0187 }
0188 }
0189 detsWithStates.clear();
0190 if (!lastTsos.isValid())
0191 return emptyContainer;
0192
0193 TrajectoryStateOnSurface secondLast = lastTsos;
0194
0195 lastTsos.rescaleError(10.0);
0196
0197 Trajectory theTraj(seed, alongMomentum);
0198
0199 navLayers.clear();
0200
0201 if (fabs(lastTsos.globalMomentum().eta()) < 1.0) {
0202
0203 navLayers = navigation()->compatibleLayers(*(lastTsos.freeState()), alongMomentum);
0204 } else if (beamhaloFlag || (theTraversingMuonFlag && theStrict1LegFlag)) {
0205
0206 navLayers = navigation()->compatibleEndcapLayers(*(lastTsos.freeState()), alongMomentum);
0207 } else {
0208 navLayers = navigation()->compatibleLayers(*(lastTsos.freeState()), alongMomentum);
0209 }
0210
0211 int DTChamberUsedBack = 0;
0212 int CSCChamberUsedBack = 0;
0213 int RPCChamberUsedBack = 0;
0214 int TotalChamberUsedBack = 0;
0215 MuonTransientTrackingRecHit::MuonRecHitContainer allUnusedHits;
0216 vector<TrajectoryMeasurement> measL;
0217
0218 LogTrace(category_) << "Begin forward fit " << navLayers.size();
0219
0220 for (vector<const DetLayer*>::const_iterator rnxtlayer = navLayers.begin(); rnxtlayer != navLayers.end();
0221 ++rnxtlayer) {
0222 LogTrace(category_) << "new layer ";
0223 measL.clear();
0224 LogTrace(category_) << debug.dumpMuonId((*rnxtlayer)->basicComponents().front()->geographicalId())
0225 << debug.dumpLayer(*rnxtlayer);
0226 LogTrace(category_) << "from lastTsos " << lastTsos.globalMomentum() << " at " << lastTsos.globalPosition();
0227
0228 measL = findBestMeasurements(*rnxtlayer, lastTsos, propagator(), (updator()->estimator()));
0229
0230 if (measL.empty() && (fabs(theService->magneticField()->inTesla(GlobalPoint(0, 0, 0)).z()) < 0.01) &&
0231 (theService->propagator("StraightLinePropagator").isValid())) {
0232 LogTrace(category_) << "try straight line propagator ";
0233 measL = findBestMeasurements(
0234 *rnxtlayer, lastTsos, &*theService->propagator("StraightLinePropagator"), (updator()->estimator()));
0235 }
0236 if (measL.empty())
0237 continue;
0238
0239 for (vector<TrajectoryMeasurement>::const_iterator theMeas = measL.begin(); theMeas != measL.end(); ++theMeas) {
0240 pair<bool, TrajectoryStateOnSurface> result = updator()->update((&*theMeas), theTraj, propagator());
0241
0242 if (result.first) {
0243 LogTrace(category_) << "update ok ";
0244 incrementChamberCounters(
0245 (*rnxtlayer), DTChamberUsedBack, CSCChamberUsedBack, RPCChamberUsedBack, TotalChamberUsedBack);
0246 secondLast = lastTsos;
0247 if ((!theTraj.empty()) && result.second.isValid()) {
0248 lastTsos = result.second;
0249 LogTrace(category_) << "get new lastTsos here " << lastTsos.globalMomentum() << " at "
0250 << lastTsos.globalPosition();
0251 } else if ((theMeas)->predictedState().isValid())
0252 lastTsos = (theMeas)->predictedState();
0253 }
0254 }
0255 }
0256 measL.clear();
0257 while (!theTraj.empty()) {
0258 theTraj.pop();
0259 }
0260
0261 if (!theTraj.isValid() || TotalChamberUsedBack < 2 || (DTChamberUsedBack + CSCChamberUsedBack) == 0 ||
0262 !lastTsos.isValid()) {
0263 return emptyContainer;
0264 }
0265
0266
0267 DTChamberUsedBack = 0;
0268 CSCChamberUsedBack = 0;
0269 RPCChamberUsedBack = 0;
0270 TotalChamberUsedBack = 0;
0271
0272 Trajectory myTraj(seed, oppositeToMomentum);
0273
0274
0275
0276 GlobalPoint lastPos = lastTsos.globalPosition();
0277 GlobalPoint secondLastPos = secondLast.globalPosition();
0278 GlobalVector momDir = secondLastPos - lastPos;
0279
0280 if (lastPos.basicVector().dot(momDir.basicVector()) > 0) {
0281
0282 theBKUpdator->setFitDirection(insideOut);
0283 } else
0284 theBKUpdator->setFitDirection(outsideIn);
0285
0286 if (fabs(lastTsos.globalMomentum().eta()) < 1.0) {
0287
0288 navLayers = navigation()->compatibleLayers(*(lastTsos.freeState()), oppositeToMomentum);
0289 } else if (beamhaloFlag || (theTraversingMuonFlag && theStrict1LegFlag)) {
0290
0291 std::reverse(navLayers.begin(), navLayers.end());
0292 } else {
0293 navLayers = navigation()->compatibleLayers(*(lastTsos.freeState()), oppositeToMomentum);
0294 }
0295
0296 LogTrace(category_) << "Begin backward refitting, with " << navLayers.size() << " layers" << endl;
0297
0298 for (vector<const DetLayer*>::const_iterator rnxtlayer = navLayers.begin(); rnxtlayer != navLayers.end();
0299 ++rnxtlayer) {
0300 measL.clear();
0301
0302 measL = findBestMeasurements(*rnxtlayer, lastTsos, propagator(), (backwardUpdator()->estimator()));
0303
0304 if (measL.empty()) {
0305 MuonTransientTrackingRecHit::MuonRecHitContainer tmpHits = theLayerMeasurements->recHits(*rnxtlayer);
0306 for (MuonRecHitContainer::const_iterator ihit = tmpHits.begin(); ihit != tmpHits.end(); ++ihit) {
0307 allUnusedHits.push_back(*ihit);
0308 }
0309 continue;
0310 }
0311
0312 for (vector<TrajectoryMeasurement>::const_iterator theMeas = measL.begin(); theMeas != measL.end(); ++theMeas) {
0313
0314 if (rnxtlayer != navLayers.begin()) {
0315 vector<const DetLayer*>::const_iterator lastlayer = rnxtlayer;
0316 lastlayer--;
0317
0318 if ((*rnxtlayer)->location() != (*lastlayer)->location()) {
0319 lastPos = lastTsos.globalPosition();
0320 GlobalPoint thisPos = (theMeas)->predictedState().globalPosition();
0321 GlobalVector momDir = thisPos - lastPos;
0322
0323
0324 if (momDir.mag() > 0.01) {
0325 if (thisPos.basicVector().dot(momDir.basicVector()) > 0) {
0326 theBKUpdator->setFitDirection(insideOut);
0327 } else
0328 theBKUpdator->setFitDirection(outsideIn);
0329 }
0330 }
0331 if (((*lastlayer)->location() == GeomDetEnumerators::endcap) &&
0332 ((*rnxtlayer)->location() == GeomDetEnumerators::endcap) &&
0333 (lastTsos.globalPosition().z() * (theMeas)->predictedState().globalPosition().z() < 0)) {
0334 theBKUpdator->setFitDirection(insideOut);
0335 }
0336 }
0337
0338
0339
0340
0341 pair<bool, TrajectoryStateOnSurface> bkresult = backwardUpdator()->update((&*theMeas), myTraj, propagator());
0342
0343 if (bkresult.first) {
0344 incrementChamberCounters(
0345 (*rnxtlayer), DTChamberUsedBack, CSCChamberUsedBack, RPCChamberUsedBack, TotalChamberUsedBack);
0346
0347 if (theTraversingMuonFlag) {
0348 MuonRecHitContainer tmpUnusedHits = unusedHits(*rnxtlayer, *theMeas);
0349 allUnusedHits.insert(allUnusedHits.end(), tmpUnusedHits.begin(), tmpUnusedHits.end());
0350 }
0351 if ((!myTraj.empty()) && bkresult.second.isValid())
0352 lastTsos = bkresult.second;
0353 else if ((theMeas)->predictedState().isValid())
0354 lastTsos = (theMeas)->predictedState();
0355 }
0356 }
0357 }
0358
0359 if ((!myTraj.isValid()) || (myTraj.empty()) || ((selfDuplicate(myTraj))) || TotalChamberUsedBack < 2 ||
0360 (DTChamberUsedBack + CSCChamberUsedBack) < 1) {
0361 return emptyContainer;
0362 }
0363
0364 if (theTraversingMuonFlag && (allUnusedHits.size() >= 2)) {
0365
0366 LogTrace(category_) << "Building trajectory in second hemisphere...";
0367 buildSecondHalf(myTraj);
0368
0369
0370 if (theStrict1LegFlag && !utilities()->isTraversing(myTraj)) {
0371 return emptyContainer;
0372 }
0373 } else if (theStrict1LegFlag && theTraversingMuonFlag) {
0374 return emptyContainer;
0375 }
0376
0377 LogTrace(category_) << " traj ok ";
0378
0379
0380 if (beamhaloFlag)
0381 estimateDirection(myTraj);
0382 if (myTraj.empty())
0383 return emptyContainer;
0384
0385
0386 vector<Trajectory> smoothed = theSmoother->trajectories(myTraj);
0387
0388 if (!smoothed.empty() && smoothed.front().foundHits() > 3) {
0389 LogTrace(category_) << " Smoothed successfully.";
0390 myTraj = smoothed.front();
0391 } else {
0392 LogTrace(category_) << " Smooth failed.";
0393 }
0394
0395 LogTrace(category_) << "first " << myTraj.firstMeasurement().updatedState() << "\n last "
0396 << myTraj.lastMeasurement().updatedState();
0397 if (myTraj.direction() == alongMomentum)
0398 LogTrace(category_) << "alongMomentum";
0399 else if (myTraj.direction() == oppositeToMomentum)
0400 LogTrace(category_) << "oppositeMomentum";
0401 else
0402 LogTrace(category_) << "anyDirection";
0403
0404 if (!beamhaloFlag) {
0405 if (myTraj.lastMeasurement().updatedState().globalMomentum().y() > 0) {
0406 LogTrace(category_) << "flip trajectory ";
0407 flipTrajectory(myTraj);
0408 }
0409
0410 if ((myTraj.direction() == alongMomentum && (myTraj.firstMeasurement().updatedState().globalPosition().y() <
0411 myTraj.lastMeasurement().updatedState().globalPosition().y())) ||
0412 (myTraj.direction() == oppositeToMomentum && (myTraj.firstMeasurement().updatedState().globalPosition().y() >
0413 myTraj.lastMeasurement().updatedState().globalPosition().y()))) {
0414 LogTrace(category_) << "reverse propagation direction";
0415 reverseTrajectoryPropagationDirection(myTraj);
0416 }
0417 }
0418
0419 if (!myTraj.isValid())
0420 return emptyContainer;
0421
0422
0423 PropagationDirection dir = myTraj.direction();
0424 GlobalVector dirFromPos = myTraj.measurements().back().recHit()->globalPosition() -
0425 myTraj.measurements().front().recHit()->globalPosition();
0426
0427 if (theStrict1LegFlag && !utilities()->isTraversing(myTraj)) {
0428 return emptyContainer;
0429 }
0430
0431 LogTrace(category_) << "last hit " << myTraj.measurements().back().recHit()->globalPosition() << endl;
0432 LogTrace(category_) << "first hit " << myTraj.measurements().front().recHit()->globalPosition() << endl;
0433
0434 LogTrace(category_) << "last tsos " << myTraj.measurements().back().updatedState().globalPosition() << " mom "
0435 << myTraj.measurements().back().updatedState().globalMomentum() << endl;
0436 LogTrace(category_) << "first tsos " << myTraj.measurements().front().updatedState().globalPosition() << " mom "
0437 << myTraj.measurements().front().updatedState().globalMomentum() << endl;
0438
0439 PropagationDirection propDir =
0440 (dirFromPos.basicVector().dot(myTraj.firstMeasurement().updatedState().globalMomentum().basicVector()) > 0)
0441 ? alongMomentum
0442 : oppositeToMomentum;
0443 LogTrace(category_) << " dir " << dir << " propDir " << propDir << endl;
0444
0445 LogTrace(category_) << "chi2 " << myTraj.chiSquared() << endl;
0446
0447 if (dir != propDir) {
0448 LogTrace(category_) << "reverse propagation direction ";
0449 reverseTrajectoryPropagationDirection(myTraj);
0450 }
0451 if (myTraj.empty())
0452 return emptyContainer;
0453
0454 navLayers.clear();
0455 TrajectoryContainer ret;
0456 ret.reserve(1);
0457 ret.emplace_back(std::make_unique<Trajectory>(myTraj));
0458 return ret;
0459 }
0460
0461
0462
0463
0464 MuonTransientTrackingRecHit::MuonRecHitContainer CosmicMuonTrajectoryBuilder::unusedHits(
0465 const DetLayer* layer, const TrajectoryMeasurement& meas) const {
0466 MuonTransientTrackingRecHit::MuonRecHitContainer tmpHits = theLayerMeasurements->recHits(layer);
0467 MuonRecHitContainer result;
0468 for (MuonRecHitContainer::const_iterator ihit = tmpHits.begin(); ihit != tmpHits.end(); ++ihit) {
0469 if ((*ihit)->geographicalId() != meas.recHit()->geographicalId()) {
0470 result.push_back(*ihit);
0471 LogTrace(category_) << "Unused hit: " << (*ihit)->globalPosition() << endl;
0472 }
0473 }
0474
0475 return result;
0476 }
0477
0478
0479
0480
0481 void CosmicMuonTrajectoryBuilder::build(const TrajectoryStateOnSurface& ts,
0482 const NavigationDirection& startingDir,
0483 Trajectory& traj) {
0484 if (!ts.isValid())
0485 return;
0486
0487 const FreeTrajectoryState* fts = ts.freeState();
0488 if (!fts)
0489 return;
0490
0491 vector<const DetLayer*> navLayers;
0492
0493 if (fabs(fts->momentum().basicVector().eta()) < 1.0) {
0494
0495 if (fts->position().basicVector().dot(fts->momentum().basicVector()) > 0) {
0496 navLayers = navigation()->compatibleLayers((*fts), alongMomentum);
0497 } else {
0498 navLayers = navigation()->compatibleLayers((*fts), oppositeToMomentum);
0499 }
0500
0501 } else if (theTraversingMuonFlag && theStrict1LegFlag) {
0502
0503 if (fts->position().basicVector().dot(fts->momentum().basicVector()) > 0) {
0504 navLayers = navigation()->compatibleEndcapLayers((*fts), alongMomentum);
0505 } else {
0506 navLayers = navigation()->compatibleEndcapLayers((*fts), oppositeToMomentum);
0507 }
0508 } else {
0509 if (fts->position().basicVector().dot(fts->momentum().basicVector()) > 0) {
0510 navLayers = navigation()->compatibleLayers((*fts), alongMomentum);
0511 } else {
0512 navLayers = navigation()->compatibleLayers((*fts), oppositeToMomentum);
0513 }
0514 }
0515
0516 if (navLayers.empty())
0517 return;
0518
0519 theBKUpdator->setFitDirection(startingDir);
0520
0521 int DTChamberUsedBack = 0;
0522 int CSCChamberUsedBack = 0;
0523 int RPCChamberUsedBack = 0;
0524 int TotalChamberUsedBack = 0;
0525
0526 TrajectoryStateOnSurface lastTsos = (traj.lastMeasurement().updatedState().globalPosition().y() <
0527 traj.firstMeasurement().updatedState().globalPosition().y())
0528 ? propagatorAlong()->propagate((*fts), navLayers.front()->surface())
0529 : propagatorOpposite()->propagate((*fts), navLayers.front()->surface());
0530
0531 if (!lastTsos.isValid()) {
0532 LogTrace(category_) << "propagation failed from fts to inner cylinder";
0533 return;
0534 }
0535 LogTrace(category_) << "tsos " << lastTsos.globalPosition();
0536 lastTsos.rescaleError(10.);
0537 vector<TrajectoryMeasurement> measL;
0538 for (vector<const DetLayer*>::const_iterator rnxtlayer = navLayers.begin(); rnxtlayer != navLayers.end();
0539 ++rnxtlayer) {
0540 measL.clear();
0541 measL = findBestMeasurements(*rnxtlayer, lastTsos, propagator(), (backwardUpdator()->estimator()));
0542
0543 if (measL.empty())
0544 continue;
0545
0546 for (vector<TrajectoryMeasurement>::const_iterator theMeas = measL.begin(); theMeas != measL.end(); ++theMeas) {
0547 pair<bool, TrajectoryStateOnSurface> bkresult = backwardUpdator()->update((&*theMeas), traj, propagator());
0548 if (bkresult.first) {
0549 LogTrace(category_) << "update ok : " << (theMeas)->recHit()->globalPosition();
0550
0551 incrementChamberCounters(
0552 (*rnxtlayer), DTChamberUsedBack, CSCChamberUsedBack, RPCChamberUsedBack, TotalChamberUsedBack);
0553
0554 if ((!traj.empty()) && bkresult.second.isValid())
0555 lastTsos = bkresult.second;
0556 else if ((theMeas)->predictedState().isValid())
0557 lastTsos = (theMeas)->predictedState();
0558 }
0559 }
0560 }
0561 navLayers.clear();
0562 updator()->makeFirstTime();
0563 backwardUpdator()->makeFirstTime();
0564
0565 measL.clear();
0566
0567 return;
0568 }
0569
0570
0571
0572
0573
0574 void CosmicMuonTrajectoryBuilder::buildSecondHalf(Trajectory& traj) {
0575 if ((traj.firstMeasurement().recHit()->globalPosition().perp() <
0576 traj.lastMeasurement().recHit()->globalPosition().perp())) {
0577 LogTrace(category_) << "inside-out: reverseTrajectory";
0578 reverseTrajectory(traj);
0579 }
0580 if (traj.empty())
0581 return;
0582 TrajectoryStateOnSurface tsos = traj.lastMeasurement().updatedState();
0583 if (!tsos.isValid())
0584 tsos = traj.lastMeasurement().predictedState();
0585 LogTrace(category_) << "last tsos on traj: pos: " << tsos.globalPosition() << " mom: " << tsos.globalMomentum();
0586 if (!tsos.isValid()) {
0587 LogTrace(category_) << "last tsos on traj invalid";
0588 return;
0589 }
0590
0591 build(intermediateState(tsos), insideOut, traj);
0592
0593 return;
0594 }
0595
0596
0597
0598
0599 TrajectoryStateOnSurface CosmicMuonTrajectoryBuilder::intermediateState(const TrajectoryStateOnSurface& tsos) const {
0600 PerpendicularBoundPlaneBuilder planeBuilder;
0601 GlobalPoint pos(0.0, 0.0, 0.0);
0602 BoundPlane* SteppingPlane = planeBuilder(pos, tsos.globalDirection());
0603
0604 TrajectoryStateOnSurface predTsos = propagator()->propagate(tsos, *SteppingPlane);
0605 if (predTsos.isValid())
0606 LogTrace(category_) << "intermediateState: a intermediate state: pos: " << predTsos.globalPosition()
0607 << "mom: " << predTsos.globalMomentum();
0608
0609 return predTsos;
0610 }
0611
0612
0613
0614
0615 void CosmicMuonTrajectoryBuilder::selectHits(MuonTransientTrackingRecHit::MuonRecHitContainer& hits) const {
0616 if (hits.size() < 2)
0617 return;
0618
0619 MuonRecHitContainer tmp;
0620 vector<bool> keep(hits.size(), true);
0621 int i(0);
0622 int j(0);
0623
0624 for (MuonRecHitContainer::const_iterator ihit = hits.begin(); ihit != hits.end(); ++ihit) {
0625 if (!keep[i]) {
0626 i++;
0627 continue;
0628 };
0629 j = i + 1;
0630 for (MuonRecHitContainer::const_iterator ihit2 = ihit + 1; ihit2 != hits.end(); ++ihit2) {
0631 if (!keep[j]) {
0632 j++;
0633 continue;
0634 }
0635 if ((*ihit)->geographicalId() == (*ihit2)->geographicalId()) {
0636 if ((*ihit)->dimension() > (*ihit2)->dimension()) {
0637 keep[j] = false;
0638 } else if ((*ihit)->dimension() < (*ihit2)->dimension()) {
0639 keep[i] = false;
0640 } else {
0641 if ((*ihit)->transientHits().size() > (*ihit2)->transientHits().size()) {
0642 keep[j] = false;
0643 } else if ((*ihit)->transientHits().size() < (*ihit2)->transientHits().size()) {
0644 keep[i] = false;
0645 } else if ((*ihit)->degreesOfFreedom() != 0 && (*ihit2)->degreesOfFreedom() != 0) {
0646 if (((*ihit)->chi2() / (*ihit)->degreesOfFreedom()) > ((*ihit2)->chi2() / (*ihit)->degreesOfFreedom()))
0647 keep[i] = false;
0648 else
0649 keep[j] = false;
0650 }
0651 }
0652 }
0653 j++;
0654 }
0655 i++;
0656 }
0657
0658 i = 0;
0659 for (MuonRecHitContainer::const_iterator ihit = hits.begin(); ihit != hits.end(); ++ihit) {
0660 if (keep[i])
0661 tmp.push_back(*ihit);
0662 i++;
0663 }
0664
0665 hits.clear();
0666 hits.swap(tmp);
0667
0668 return;
0669 }
0670
0671
0672
0673
0674 bool CosmicMuonTrajectoryBuilder::selfDuplicate(const Trajectory& traj) const {
0675 TransientTrackingRecHit::ConstRecHitContainer const& hits = traj.recHits();
0676
0677 if (traj.empty())
0678 return true;
0679
0680 bool result = false;
0681 for (ConstRecHitContainer::const_iterator ir = hits.begin(); ir != hits.end(); ir++) {
0682 if (!(*ir)->isValid())
0683 continue;
0684 for (ConstRecHitContainer::const_iterator ir2 = ir + 1; ir2 != hits.end(); ir2++) {
0685 if (!(*ir2)->isValid())
0686 continue;
0687 if ((*ir) == (*ir2))
0688 result = true;
0689 }
0690 }
0691
0692 return result;
0693 }
0694
0695
0696
0697
0698
0699 void CosmicMuonTrajectoryBuilder::reverseTrajectory(Trajectory& traj) const {
0700 PropagationDirection newDir =
0701 (traj.firstMeasurement().recHit()->globalPosition().y() < traj.lastMeasurement().recHit()->globalPosition().y())
0702 ? oppositeToMomentum
0703 : alongMomentum;
0704 Trajectory newTraj(traj.seed(), newDir);
0705
0706
0707
0708
0709
0710
0711
0712
0713
0714 std::vector<TrajectoryMeasurement> const& meas = traj.measurements();
0715 for (auto itm = meas.rbegin(); itm != meas.rend(); ++itm) {
0716 newTraj.push(*itm);
0717 }
0718 traj = newTraj;
0719 }
0720
0721
0722
0723
0724 void CosmicMuonTrajectoryBuilder::flipTrajectory(Trajectory& traj) const {
0725 TrajectoryStateOnSurface lastTSOS = traj.lastMeasurement().updatedState();
0726 if (!lastTSOS.isValid()) {
0727 LogTrace(category_) << "Error: last TrajectoryState invalid.";
0728 }
0729 TransientTrackingRecHit::ConstRecHitContainer hits = traj.recHits();
0730 std::reverse(hits.begin(), hits.end());
0731
0732 LogTrace(category_) << "last tsos before flipping " << lastTSOS;
0733 utilities()->reverseDirection(lastTSOS, &*theService->magneticField());
0734 LogTrace(category_) << "last tsos after flipping " << lastTSOS;
0735
0736 vector<Trajectory> refittedback = theSmoother->fit(traj.seed(), hits, lastTSOS);
0737 if (refittedback.empty()) {
0738 LogTrace(category_) << "flipTrajectory fail. " << endl;
0739 return;
0740 }
0741 LogTrace(category_) << "flipTrajectory: first " << refittedback.front().firstMeasurement().updatedState()
0742 << "\nflipTrajectory: last " << refittedback.front().lastMeasurement().updatedState();
0743
0744 traj = refittedback.front();
0745
0746 return;
0747 }
0748
0749
0750
0751
0752 void CosmicMuonTrajectoryBuilder::reverseTrajectoryPropagationDirection(Trajectory& traj) const {
0753 if (traj.direction() == anyDirection)
0754 return;
0755 PropagationDirection newDir = (traj.direction() == alongMomentum) ? oppositeToMomentum : alongMomentum;
0756 Trajectory newTraj(traj.seed(), newDir);
0757 const std::vector<TrajectoryMeasurement>& meas = traj.measurements();
0758
0759 for (std::vector<TrajectoryMeasurement>::const_iterator itm = meas.begin(); itm != meas.end(); ++itm) {
0760 newTraj.push(*itm);
0761 }
0762
0763 while (!traj.empty()) {
0764 traj.pop();
0765 }
0766
0767 traj = newTraj;
0768 }
0769
0770
0771
0772
0773 void CosmicMuonTrajectoryBuilder::estimateDirection(Trajectory& traj) const {
0774 TransientTrackingRecHit::ConstRecHitContainer hits = traj.recHits();
0775
0776 TrajectoryStateOnSurface firstTSOS = traj.firstMeasurement().updatedState();
0777
0778 TrajectoryStateOnSurface lastTSOS = traj.lastMeasurement().updatedState();
0779
0780 if (!firstTSOS.isValid() || !lastTSOS.isValid())
0781 return;
0782
0783 LogTrace(category_) << "Two ends of the traj " << firstTSOS.globalPosition() << ", " << lastTSOS.globalPosition();
0784
0785 LogTrace(category_) << "Their mom: " << firstTSOS.globalMomentum() << ", " << lastTSOS.globalMomentum();
0786
0787 LogTrace(category_) << "Their mom eta: " << firstTSOS.globalMomentum().eta() << ", "
0788 << lastTSOS.globalMomentum().eta();
0789
0790
0791
0792
0793 if (fabs(firstTSOS.globalMomentum().eta()) > fabs(lastTSOS.globalMomentum().eta())) {
0794 vector<Trajectory> refitted = theSmoother->trajectories(traj.seed(), hits, firstTSOS);
0795 if (!refitted.empty())
0796 traj = refitted.front();
0797
0798 } else {
0799 std::reverse(hits.begin(), hits.end());
0800 utilities()->reverseDirection(lastTSOS, &*theService->magneticField());
0801 vector<Trajectory> refittedback = theSmoother->trajectories(traj.seed(), hits, lastTSOS);
0802 if (!refittedback.empty())
0803 traj = refittedback.front();
0804 }
0805
0806 return;
0807 }
0808
0809
0810
0811
0812 void CosmicMuonTrajectoryBuilder::getDirectionByTime(Trajectory& traj) const {
0813 TransientTrackingRecHit::ConstRecHitContainer hits = traj.recHits();
0814 LogTrace(category_) << "getDirectionByTime" << endl;
0815 for (TransientTrackingRecHit::ConstRecHitContainer::const_iterator ir = hits.begin(); ir != hits.end(); ir++) {
0816 if (!(*ir)->isValid()) {
0817 LogTrace(category_) << "invalid RecHit" << endl;
0818 continue;
0819 }
0820
0821 const GlobalPoint& pos = (*ir)->globalPosition();
0822 LogTrace(category_) << "pos" << pos << "radius " << pos.perp() << " dim " << (*ir)->dimension() << " det "
0823 << (*ir)->det()->geographicalId().det() << " sub det " << (*ir)->det()->subDetector() << endl;
0824
0825 if ((*ir)->det()->geographicalId().det() == 2 && (*ir)->det()->subDetector() == 6) {
0826
0827
0828 CSCRecHit2DCollection::range thisrange = cschits_->get(CSCDetId((*ir)->geographicalId()));
0829 for (CSCRecHit2DCollection::const_iterator rechit = thisrange.first; rechit != thisrange.second; ++rechit) {
0830 if ((*rechit).isValid())
0831 LogTrace(category_) << "csc from collection tpeak " << (*rechit).tpeak();
0832 }
0833 }
0834 if ((*ir)->det()->geographicalId().det() == 2 && (*ir)->det()->subDetector() == 7) {
0835
0836
0837 DTRecHitCollection::range thisrange = dthits_->get(DTLayerId((*ir)->geographicalId()));
0838 for (DTRecHitCollection::const_iterator rechit = thisrange.first; rechit != thisrange.second; ++rechit) {
0839 if ((*rechit).isValid())
0840 LogTrace(category_) << "dt from collection digitime " << (*rechit).digiTime();
0841 }
0842 }
0843 }
0844
0845 return;
0846 }
0847
0848
0849
0850
0851 std::vector<TrajectoryMeasurement> CosmicMuonTrajectoryBuilder::findBestMeasurements(
0852 const DetLayer* layer,
0853 const TrajectoryStateOnSurface& tsos,
0854 const Propagator* propagator,
0855 const MeasurementEstimator* estimator) {
0856 std::vector<TrajectoryMeasurement> result;
0857 std::vector<TrajectoryMeasurement> measurements;
0858
0859 if (layer->hasGroups()) {
0860 std::vector<TrajectoryMeasurementGroup> measurementGroups =
0861 theLayerMeasurements->groupedMeasurements(layer, tsos, *propagator, *estimator);
0862
0863 for (std::vector<TrajectoryMeasurementGroup>::const_iterator tmGroupItr = measurementGroups.begin();
0864 tmGroupItr != measurementGroups.end();
0865 ++tmGroupItr) {
0866 measurements = tmGroupItr->measurements();
0867 const TrajectoryMeasurement* bestMeasurement =
0868 theBestMeasurementFinder->findBestMeasurement(measurements, propagator);
0869
0870 if (bestMeasurement)
0871 result.push_back(*bestMeasurement);
0872 }
0873 } else {
0874 measurements = theLayerMeasurements->measurements(layer, tsos, *propagator, *estimator);
0875 const TrajectoryMeasurement* bestMeasurement =
0876 theBestMeasurementFinder->findBestMeasurement(measurements, propagator);
0877
0878 if (bestMeasurement)
0879 result.push_back(*bestMeasurement);
0880 }
0881 measurements.clear();
0882
0883 return result;
0884 }
0885
0886
0887
0888
0889 void CosmicMuonTrajectoryBuilder::incrementChamberCounters(
0890 const DetLayer* layer, int& dtChambers, int& cscChambers, int& rpcChambers, int& totalChambers) {
0891 if (layer->subDetector() == GeomDetEnumerators::DT)
0892 dtChambers++;
0893 else if (layer->subDetector() == GeomDetEnumerators::CSC)
0894 cscChambers++;
0895 else if (layer->subDetector() == GeomDetEnumerators::RPCBarrel ||
0896 layer->subDetector() == GeomDetEnumerators::RPCEndcap)
0897 rpcChambers++;
0898 totalChambers++;
0899 }
0900
0901
0902
0903
0904 double CosmicMuonTrajectoryBuilder::t0(const DTRecSegment4D* dtseg) const {
0905 if ((dtseg == nullptr) || (!dtseg->hasPhi()))
0906 return 0;
0907
0908 double result = 0;
0909 if (dtseg->phiSegment() == nullptr)
0910 return 0;
0911 int phiHits = dtseg->phiSegment()->specificRecHits().size();
0912 LogTrace(category_) << "phiHits " << phiHits;
0913 if (phiHits > 5) {
0914 if (dtseg->phiSegment()->ist0Valid())
0915 result = dtseg->phiSegment()->t0();
0916 if (dtseg->phiSegment()->ist0Valid()) {
0917 LogTrace(category_) << " Phi t0: " << dtseg->phiSegment()->t0() << " hits: " << phiHits;
0918 } else {
0919 LogTrace(category_) << " Phi t0 is invalid: " << dtseg->phiSegment()->t0() << " hits: " << phiHits;
0920 }
0921 }
0922
0923 return result;
0924 }
0925
0926
0927
0928
0929 PropagationDirection CosmicMuonTrajectoryBuilder::checkDirectionByT0(const DTRecSegment4D* dtseg1,
0930 const DTRecSegment4D* dtseg2) const {
0931 LogTrace(category_) << "comparing dtseg: " << dtseg1 << " " << dtseg2 << endl;
0932 if (dtseg1 == dtseg2 || t0(dtseg1) == t0(dtseg2))
0933 return anyDirection;
0934
0935 PropagationDirection result = (t0(dtseg1) < t0(dtseg2)) ? alongMomentum : oppositeToMomentum;
0936
0937 return result;
0938 }