Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-10-25 10:01:08 

0001 #include "RecoMuon/CosmicMuonProducer/interface/CosmicMuonTrajectoryBuilder.h"
0002 /** \file CosmicMuonTrajectoryBuilder
0003  *
0004  *  class to build trajectories of cosmic muons and beam-halo muons
0005  *
0006  *
0007  *  \author Chang Liu  - Purdue Univeristy
0008  */
0009 
0010 #include "FWCore/Framework/interface/Frameworkfwd.h"
0011 #include "FWCore/Framework/interface/MakerMacros.h"
0012 
0013 /* Collaborating Class Header */
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   //  if(enableDTMeasurement)
0054   InputTag DTRecSegmentLabel = par.getParameter<InputTag>("DTRecSegmentLabel");
0055 
0056   //  if(enableCSCMeasurement)
0057   InputTag CSCRecSegmentLabel = par.getParameter<InputTag>("CSCRecSegmentLabel");
0058 
0059   //  if(enableRPCMeasurement)
0060   InputTag RPCRecSegmentLabel = par.getParameter<InputTag>("RPCRecSegmentLabel");
0061 
0062   //  if(enableGEMMeasurement)
0063   //  InputTag GEMRecSegmentLabel = par.getParameter<InputTag>("GEMRecSegmentLabel");
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;  // new DirectMuonNavigation(theService->detLayerGeometry());
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   // DetLayer Geometry
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     //DT
0155     navLayers = navigation()->compatibleLayers(*(lastTsos.freeState()), alongMomentum);
0156   } else if (beamhaloFlag || (theTraversingMuonFlag && theStrict1LegFlag)) {
0157     //CSC
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     // get the updated TSOS
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     //DT
0203     navLayers = navigation()->compatibleLayers(*(lastTsos.freeState()), alongMomentum);
0204   } else if (beamhaloFlag || (theTraversingMuonFlag && theStrict1LegFlag)) {
0205     //CSC
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   // if got good trajectory, then do backward refitting
0267   DTChamberUsedBack = 0;
0268   CSCChamberUsedBack = 0;
0269   RPCChamberUsedBack = 0;
0270   TotalChamberUsedBack = 0;
0271 
0272   Trajectory myTraj(seed, oppositeToMomentum);
0273 
0274   // set starting navigation direction for MuonTrajectoryUpdator
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     //      LogTrace("CosmicMuonTrajectoryBuilder")<<"Fit direction changed to insideOut";
0282     theBKUpdator->setFitDirection(insideOut);
0283   } else
0284     theBKUpdator->setFitDirection(outsideIn);
0285 
0286   if (fabs(lastTsos.globalMomentum().eta()) < 1.0) {
0287     //DT
0288     navLayers = navigation()->compatibleLayers(*(lastTsos.freeState()), oppositeToMomentum);
0289   } else if (beamhaloFlag || (theTraversingMuonFlag && theStrict1LegFlag)) {
0290     //CSC
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       // if the part change, we need to reconsider the fit direction
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           //         LogTrace("CosmicMuonTrajectoryBuilder")<<"momDir "<<momDir;
0323 
0324           if (momDir.mag() > 0.01) {  //if lastTsos is on the surface, no need
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       //       if (theBKUpdator->fitDirection() == insideOut)
0339       //          LogTrace("CosmicMuonTrajectoryBuilder")<<"Fit direction insideOut";
0340       //       else LogTrace("CosmicMuonTrajectoryBuilder")<<"Fit direction outsideIn";
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     //      LogTrace(category_)<<utilities()->print(allUnusedHits);
0366     LogTrace(category_) << "Building trajectory in second hemisphere...";
0367     buildSecondHalf(myTraj);
0368     // check if traversing trajectory has hits in both hemispheres
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   //     getDirectionByTime(myTraj);
0380   if (beamhaloFlag)
0381     estimateDirection(myTraj);
0382   if (myTraj.empty())
0383     return emptyContainer;
0384 
0385   // try to smooth it
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   //  getDirectionByTime(myTraj);
0419   if (!myTraj.isValid())
0420     return emptyContainer;
0421 
0422   // check direction agree with position!
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 // continue to build a trajectory starting from given trajectory state
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     //DT
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     //CSC
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 // build trajectory in second hemisphere with pattern
0572 // recognition starting from an intermediate state
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       }  // if same geomid
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 // reverse a trajectory without refitting
0697 // this can only be used for cosmic muons that come from above
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   /* does not work in gcc4.8?)
0707   std::vector<TrajectoryMeasurement> & meas = traj.measurements();
0708   for (auto itm = meas.rbegin(); itm != meas.rend(); ++itm ) {
0709     newTraj.push(std::move(*itm));
0710   }
0711   traj = std::move(newTraj);
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 // reverse a trajectory momentum direction and then refit
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 // guess the direction by normalized chi2
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   // momentum eta can be used to estimate direction
0791   // the beam-halo muon seems enter with a larger |eta|
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 // get direction from timing information of rechits and segments
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       //      const CSCRecHit2D* iCSC = dynamic_cast<const CSCRecHit2D*>(&**ir);
0827       //      if (iCSC) LogTrace(category_)<<"csc from cast tpeak "<<iCSC->tpeak();
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       //      const DTRecHit1D* iDT = dynamic_cast<const DTRecHit1D*>(&**ir);
0836       //      if (iDT) LogTrace(category_)<<"dt digitime "<<iDT->digiTime();
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   // timing information
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 }
[ 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!