CosmicMuonProducer/src/GlobalCosmicMuonTrajectoryBuilder.cc

0001 /**
0002  *  Class: GlobalCosmicMuonTrajectoryBuilder
0003  *
0004  *  \author Chang Liu  -  Purdue University <Chang.Liu@cern.ch>
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 // constructor
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 // destructor
0050 //
0051
0052 GlobalCosmicMuonTrajectoryBuilder::~GlobalCosmicMuonTrajectoryBuilder() {
0053   if (theSmoother)
0054     delete theSmoother;
0055   if (theTrackMatcher)
0056     delete theTrackMatcher;
0057 }
0058
0059 //
0060 // set Event
0061 //
0062 void GlobalCosmicMuonTrajectoryBuilder::setEvent(const edm::Event& event) {
0063   event.getByToken(theTkTrackToken, theTrackerTracks);
0064
0065   //  edm::Handle<std::vector<Trajectory> > handleTrackerTrajs;
0066   //  if ( event.getByLabel(theTkTrackLabel,handleTrackerTrajs) && handleTrackerTrajs.isValid() ) {
0067   //      tkTrajsAvailable = true;
0068   //      allTrackerTrajs = &*handleTrackerTrajs;
0069   //      LogInfo("GlobalCosmicMuonTrajectoryBuilder")
0070   //    << "Tk Trajectories Found! " << endl;
0071   //  } else {
0072   //      LogInfo("GlobalCosmicMuonTrajectoryBuilder")
0073   //    << "No Tk Trajectories Found! " << endl;
0074   //      tkTrajsAvailable = false;
0075   //  }
0076
0077   theTrackerRecHitBuilder = theService->eventSetup().getHandle(theTrackerRecHitBuilderToken);
0078   theMuonRecHitBuilder = theService->eventSetup().getHandle(theMuonRecHitBuilderToken);
0079 }
0080
0081 //
0082 // reconstruct trajectories
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   //  if ( !tkTrajsAvailable ) {
0121   //     tkRecHits = getTransientRecHits(*tkTrack);
0122   //  } else {
0123   //     tkRecHits = allTrackerTrajs->front().recHits();
0124   //  }
0125
0126   ConstRecHitContainer hits;  //= tkRecHits;
0127   LogTrace(category_) << "tk RecHits: " << tkRecHits.size();
0128
0129   //  hits.insert(hits.end(), muRecHits.begin(), muRecHits.end());
0130   //  stable_sort(hits.begin(), hits.end(), DecreasingGlobalY());
0131
0132   sortHits(hits, muRecHits, tkRecHits);
0133
0134   //  LogTrace(category_)<< "Used RecHits after sort: "<<hits.size()<<endl;;
0135   //  LogTrace(category_) <<utilities()->print(hits)<<endl;
0136   //  LogTrace(category_) << "== End of Used RecHits == "<<endl;
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   // begin refitting
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);  //FIXME
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   //sort hits going from higher to lower positions
0255   if (frontTkPos.y() < backTkPos.y()) {  //check if tk hits order same direction
0256     reverse(tkHits.begin(), tkHits.end());
0257   }
0258
0259   if (frontMuPos.y() < backMuPos.y()) {
0260     reverse(muonHits.begin(), muonHits.end());
0261   }
0262
0263   //  LogTrace(category_)<< "tkHits after sort: "<<tkHits.size()<<endl;;
0264   //  LogTrace(category_) <<utilities()->print(tkHits)<<endl;
0265   //  LogTrace(category_) << "== End of tkHits == "<<endl;
0266
0267   //  LogTrace(category_)<< "muonHits after sort: "<<muonHits.size()<<endl;;
0268   //  LogTrace(category_) <<utilities()->print(muonHits)<<endl;
0269   //  LogTrace(category_)<< "== End of muonHits == "<<endl;
0270
0271   //separate muon hits into 2 different hemisphere
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   //insert track hits in correct order
0292   if (insertInMiddle) {  //if tk hits should be sandwich
0293     GlobalPoint jointpointpos = (*middlepoint)->globalPosition();
0294     LogTrace(category_) << "jointpoint " << jointpointpos << endl;
0295     if ((frontTkPos - jointpointpos).mag() >
0296         (backTkPos - jointpointpos).mag()) {  //check if tk hits order same direction
0297       reverse(tkHits.begin(), tkHits.end());
0298     }
0299     muonHits.insert(middlepoint + 1, tkHits.begin(), tkHits.end());
0300     hits = muonHits;
0301   } else {                                  // append at one end
0302     if (frontTkPos.y() < frontMuPos.y()) {  //insert at the end
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 {  //insert at the beginning
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   //build tracker TrackCands and pick the best match if size greater than 2
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   // now if only 1 tracker tracks, return it
0364   if (tkTrackCands.size() <= 1) {
0365     return tkTrackCands;
0366   }
0367
0368   // if there're many tracker tracks
0369
0370   // if muon is only on one side
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     // if there're tracker tracks totally on the other half
0381     // and there're tracker tracks on the same half
0382     // remove the tracks on the other half
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       //if all tk tracks on the other side, still keep them
0429       result = tkTrackCands;
0430     }
0431   } else {  // muon is traversing
0432     result = tkTrackCands;
0433   }
0434
0435   // match muCand to tkTrackCands
0436   vector<TrackCand> matched_trackerTracks = theTrackMatcher->match(mu, result);
0437
0438   LogTrace(category_) << "TrackMatcher found " << matched_trackerTracks.size() << "tracker tracks matched";
0439
0440   //now pick the best matched one
0441   if (matched_trackerTracks.size() < 2) {
0442     return matched_trackerTracks;
0443   } else {
0444     // in case of more than 1 tkTrack,
0445     // select the best-one based on distance (matchOption==1)
0446     // at innermost Mu hit surface. (surfaceOption == 0)
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 }