Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:26:51 

0001 /** \file CosmicMuonSmoother
0002  *
0003  *  The class refit muon trajectories based on MuonTrackReFitter
0004  *  But 
0005  *  (1) check direction to make sure it's downward for cosmic
0006  *  (2) propagate by virtual intermediate planes if failed to ensure propagation
0007  *      within cylinders
0008  *
0009  *
0010  *  \author Chang Liu  -  Purdue University
0011  */
0012 
0013 #include "RecoMuon/CosmicMuonProducer/interface/CosmicMuonSmoother.h"
0014 
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
0017 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0018 #include "TrackingTools/TrackFitters/interface/TrajectoryStateWithArbitraryError.h"
0019 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
0020 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0021 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0022 #include "DataFormats/GeometrySurface/interface/PlaneBuilder.h"
0023 
0024 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0025 
0026 using namespace edm;
0027 using namespace std;
0028 
0029 //
0030 // constructor
0031 //
0032 CosmicMuonSmoother::CosmicMuonSmoother(const ParameterSet& par, const MuonServiceProxy* service) : theService(service) {
0033   theUpdator = new KFUpdator;
0034   theUtilities = new CosmicMuonUtilities;
0035   theEstimator = new Chi2MeasurementEstimator(200.0);
0036   thePropagatorAlongName = par.getParameter<string>("PropagatorAlong");
0037   thePropagatorOppositeName = par.getParameter<string>("PropagatorOpposite");
0038   theErrorRescaling = par.getParameter<double>("RescalingFactor");
0039 
0040   category_ = "Muon|RecoMuon|CosmicMuon|CosmicMuonSmoother";
0041 }
0042 
0043 //
0044 // destructor
0045 //
0046 CosmicMuonSmoother::~CosmicMuonSmoother() {
0047   if (theUpdator)
0048     delete theUpdator;
0049   if (theUtilities)
0050     delete theUtilities;
0051   if (theEstimator)
0052     delete theEstimator;
0053 }
0054 
0055 //
0056 // fit and smooth trajectory
0057 //
0058 Trajectory CosmicMuonSmoother::trajectory(const Trajectory& t) const {
0059   std::vector<Trajectory>&& fitted = fit(t);
0060   if (fitted.empty())
0061     return Trajectory();
0062   std::vector<Trajectory>&& smoothed = smooth(fitted);
0063   return smoothed.empty() ? Trajectory() : smoothed.front();
0064 }
0065 
0066 //
0067 // fit and smooth trajectory
0068 //
0069 std::vector<Trajectory> CosmicMuonSmoother::trajectories(const TrajectorySeed& seed,
0070                                                          const ConstRecHitContainer& hits,
0071                                                          const TrajectoryStateOnSurface& firstPredTsos) const {
0072   if (hits.empty() || !firstPredTsos.isValid())
0073     return vector<Trajectory>();
0074 
0075   LogTrace(category_) << "trajectory begin (seed hits tsos)";
0076 
0077   TrajectoryStateOnSurface firstTsos = firstPredTsos;
0078   firstTsos.rescaleError(theErrorRescaling);
0079 
0080   LogTrace(category_) << "first TSOS: " << firstTsos;
0081 
0082   vector<Trajectory> fitted = fit(seed, hits, firstTsos);
0083   LogTrace(category_) << "fitted: " << fitted.size();
0084   vector<Trajectory> smoothed = smooth(fitted);
0085   LogTrace(category_) << "smoothed: " << smoothed.size();
0086 
0087   return smoothed;
0088 }
0089 
0090 //
0091 // fit trajectory
0092 //
0093 vector<Trajectory> CosmicMuonSmoother::fit(const Trajectory& t) const {
0094   if (t.empty())
0095     return vector<Trajectory>();
0096 
0097   LogTrace(category_) << "fit begin (trajectory) ";
0098 
0099   TrajectoryStateOnSurface firstTsos = initialState(t);
0100   if (!firstTsos.isValid()) {
0101     LogTrace(category_) << "Error: firstTsos invalid. ";
0102     return vector<Trajectory>();
0103   }
0104 
0105   LogTrace(category_) << "firstTsos: " << firstTsos;
0106 
0107   ConstRecHitContainer hits = t.recHits();
0108   LogTrace(category_) << "hits: " << hits.size();
0109   LogTrace(category_) << "hit front" << hits.front()->globalPosition() << " hit back" << hits.back()->globalPosition();
0110 
0111   sortHitsAlongMom(hits, firstTsos);
0112 
0113   LogTrace(category_) << "after sorting hit front" << hits.front()->globalPosition() << " hit back"
0114                       << hits.back()->globalPosition();
0115 
0116   return fit(t.seed(), hits, firstTsos);
0117 }
0118 
0119 //
0120 // fit trajectory
0121 //
0122 vector<Trajectory> CosmicMuonSmoother::fit(const TrajectorySeed& seed,
0123                                            const ConstRecHitContainer& hits,
0124                                            const TrajectoryStateOnSurface& firstPredTsos) const {
0125   LogTrace(category_) << "fit begin (seed, hit, tsos).";
0126 
0127   if (hits.empty()) {
0128     LogTrace(category_) << "Error: empty hits container.";
0129     return vector<Trajectory>();
0130   }
0131 
0132   Trajectory myTraj(seed, alongMomentum);
0133 
0134   TrajectoryStateOnSurface predTsos(firstPredTsos);
0135   LogTrace(category_) << "first pred TSOS: " << predTsos;
0136 
0137   if (!predTsos.isValid()) {
0138     LogTrace(category_) << "Error: firstTsos invalid.";
0139     return vector<Trajectory>();
0140   }
0141   TrajectoryStateOnSurface currTsos;
0142 
0143   if (hits.front()->isValid()) {
0144     // FIXME  FIXME  CLONE !!!
0145     //  TrackingRecHit::RecHitPointer preciseHit = hits.front()->clone(predTsos);
0146     const auto& preciseHit = hits.front();
0147 
0148     LogTrace(category_) << "first hit is at det " << hits.front()->det()->surface().position();
0149 
0150     currTsos = theUpdator->update(predTsos, *preciseHit);
0151     if (!currTsos.isValid()) {
0152       edm::LogWarning(category_) << "an updated state is not valid. killing the trajectory.";
0153       return vector<Trajectory>();
0154     }
0155     myTraj.push(TrajectoryMeasurement(
0156         predTsos, currTsos, hits.front(), theEstimator->estimate(predTsos, *hits.front()).second));
0157     if (currTsos.isValid())
0158       LogTrace(category_) << "first curr TSOS: " << currTsos;
0159 
0160   } else {
0161     currTsos = predTsos;
0162     myTraj.push(TrajectoryMeasurement(predTsos, hits.front()));
0163   }
0164   //const TransientTrackingRecHit& firsthit = *hits.front();
0165 
0166   for (ConstRecHitContainer::const_iterator ihit = hits.begin() + 1; ihit != hits.end(); ++ihit) {
0167     if (!(**ihit).isValid()) {
0168       LogTrace(category_) << "Error: invalid hit.";
0169       continue;
0170     }
0171     if (currTsos.isValid() && currTsos.globalMomentum().mag2() > 1e-18f) {
0172       LogTrace(category_) << "current pos " << currTsos.globalPosition() << "mom " << currTsos.globalMomentum();
0173     } else {
0174       LogTrace(category_) << "current state invalid or momentum is too low";
0175       //logically, there's no way out: can't expect a valid result out of an invalid state
0176       LogTrace(category_) << "Input state is not valid. This loop over hits is doomed: breaking out";
0177       break;
0178     }
0179 
0180     predTsos = propagatorAlong()->propagate(currTsos, (**ihit).det()->surface());
0181     LogTrace(category_) << "predicted state propagate directly " << predTsos.isValid();
0182 
0183     if (!predTsos.isValid()) {
0184       LogTrace(category_) << "step-propagating from " << currTsos.globalPosition()
0185                           << " to position: " << (*ihit)->globalPosition();
0186       predTsos = theUtilities->stepPropagate(currTsos, (*ihit), *propagatorAlong());
0187     }
0188     if (!predTsos.isValid() && (fabs(theService->magneticField()->inTesla(GlobalPoint(0, 0, 0)).z()) < 0.01) &&
0189         (theService->propagator("StraightLinePropagator").isValid())) {
0190       LogTrace(category_) << "straight-line propagating from " << currTsos.globalPosition()
0191                           << " to position: " << (*ihit)->globalPosition();
0192       predTsos = theService->propagator("StraightLinePropagator")->propagate(currTsos, (**ihit).det()->surface());
0193     }
0194     if (predTsos.isValid()) {
0195       LogTrace(category_) << "predicted pos " << predTsos.globalPosition() << "mom " << predTsos.globalMomentum();
0196     } else {
0197       LogTrace(category_) << "predicted state invalid";
0198     }
0199     if (!predTsos.isValid() || predTsos.globalPosition().mag2() > 1.e12) {
0200       LogTrace(category_) << "Error: predTsos is still invalid or too far in forward fit.";
0201       //      return vector<Trajectory>();
0202       continue;
0203     } else if ((**ihit).isValid()) {
0204       // FIXME  FIXME  CLONE !!!
0205       // update  (FIXME!)
0206       // TransientTrackingRecHit::RecHitPointer preciseHit = (**ihit).clone(predTsos);
0207       const auto& preciseHit = *ihit;
0208 
0209       if (!preciseHit->isValid()) {
0210         currTsos = predTsos;
0211         myTraj.push(TrajectoryMeasurement(predTsos, *ihit));
0212       } else {
0213         currTsos = theUpdator->update(predTsos, *preciseHit);
0214         if (!currTsos.isValid()) {
0215           edm::LogWarning(category_) << "an updated state is not valid. killing the trajectory.";
0216           return vector<Trajectory>();
0217         }
0218         myTraj.push(TrajectoryMeasurement(
0219             predTsos, currTsos, preciseHit, theEstimator->estimate(predTsos, *preciseHit).second));
0220       }
0221     } else {
0222       currTsos = predTsos;
0223       myTraj.push(TrajectoryMeasurement(predTsos, *ihit));
0224     }
0225   }
0226 
0227   std::vector<TrajectoryMeasurement> mytms = myTraj.measurements();
0228   LogTrace(category_) << "fit result " << mytms.size();
0229   for (std::vector<TrajectoryMeasurement>::const_iterator itm = mytms.begin(); itm != mytms.end(); ++itm) {
0230     LogTrace(category_) << "updated pos " << itm->updatedState().globalPosition() << "mom "
0231                         << itm->updatedState().globalMomentum();
0232   }
0233 
0234   return vector<Trajectory>(1, myTraj);
0235 }
0236 
0237 //
0238 // smooth trajectory
0239 //
0240 vector<Trajectory> CosmicMuonSmoother::smooth(const vector<Trajectory>& tc) const {
0241   vector<Trajectory> result;
0242 
0243   for (vector<Trajectory>::const_iterator it = tc.begin(); it != tc.end(); ++it) {
0244     vector<Trajectory> smoothed = smooth(*it);
0245     result.insert(result.end(), smoothed.begin(), smoothed.end());
0246   }
0247 
0248   return result;
0249 }
0250 
0251 //
0252 // smooth trajectory
0253 //
0254 vector<Trajectory> CosmicMuonSmoother::smooth(const Trajectory& t) const {
0255   if (t.empty()) {
0256     LogTrace(category_) << "Error: smooth: empty trajectory.";
0257     return vector<Trajectory>();
0258   }
0259 
0260   Trajectory myTraj(t.seed(), oppositeToMomentum);
0261 
0262   vector<TrajectoryMeasurement> avtm = t.measurements();
0263 
0264   if (avtm.size() < 2) {
0265     LogTrace(category_) << "Error: smooth: too little TM. ";
0266     return vector<Trajectory>();
0267   }
0268 
0269   TrajectoryStateOnSurface predTsos = avtm.back().forwardPredictedState();
0270   predTsos.rescaleError(theErrorRescaling);
0271 
0272   if (!predTsos.isValid()) {
0273     LogTrace(category_) << "Error: smooth: first TSOS from back invalid. ";
0274     return vector<Trajectory>();
0275   }
0276 
0277   TrajectoryStateOnSurface currTsos;
0278 
0279   // first smoothed TrajectoryMeasurement is last fitted
0280   if (avtm.back().recHit()->isValid()) {
0281     currTsos = theUpdator->update(predTsos, (*avtm.back().recHit()));
0282     if (!currTsos.isValid()) {
0283       edm::LogWarning(category_) << "an updated state is not valid. killing the trajectory.";
0284       return vector<Trajectory>();
0285     }
0286     myTraj.push(TrajectoryMeasurement(avtm.back().forwardPredictedState(),
0287                                       predTsos,
0288                                       avtm.back().updatedState(),
0289                                       avtm.back().recHit(),
0290                                       avtm.back().estimate()  //,
0291                                       /*avtm.back().layer()*/),
0292                 avtm.back().estimate());
0293 
0294   } else {
0295     currTsos = predTsos;
0296     myTraj.push(TrajectoryMeasurement(avtm.back().forwardPredictedState(),
0297                                       avtm.back().recHit()  //,
0298                                       /*avtm.back().layer()*/));
0299   }
0300 
0301   TrajectoryStateCombiner combiner;
0302 
0303   for (vector<TrajectoryMeasurement>::reverse_iterator itm = avtm.rbegin() + 1; itm != avtm.rend() - 1; ++itm) {
0304     if (currTsos.isValid()) {
0305       LogTrace(category_) << "current pos " << currTsos.globalPosition() << "mom " << currTsos.globalMomentum();
0306     } else {
0307       LogTrace(category_) << "current state invalid";
0308     }
0309 
0310     predTsos = propagatorOpposite()->propagate(currTsos, (*itm).recHit()->det()->surface());
0311 
0312     if (!predTsos.isValid()) {
0313       LogTrace(category_) << "step-propagating from " << currTsos.globalPosition()
0314                           << " to position: " << (*itm).recHit()->globalPosition();
0315       predTsos = theUtilities->stepPropagate(currTsos, (*itm).recHit(), *propagatorOpposite());
0316     }
0317     if (predTsos.isValid()) {
0318       LogTrace(category_) << "predicted pos " << predTsos.globalPosition() << "mom " << predTsos.globalMomentum();
0319     } else {
0320       LogTrace(category_) << "predicted state invalid";
0321     }
0322 
0323     if (!predTsos.isValid()) {
0324       LogTrace(category_) << "Error: predTsos is still invalid backward smooth.";
0325       return vector<Trajectory>();
0326       //  continue;
0327     } else if ((*itm).recHit()->isValid()) {
0328       //update
0329       currTsos = theUpdator->update(predTsos, (*(*itm).recHit()));
0330       if (!currTsos.isValid()) {
0331         edm::LogWarning(category_) << "an updated state is not valid. killing the trajectory.";
0332         return vector<Trajectory>();
0333       }
0334       TrajectoryStateOnSurface combTsos = combiner(predTsos, (*itm).forwardPredictedState());
0335       if (!combTsos.isValid()) {
0336         LogTrace(category_) << "Error: smooth: combining pred TSOS failed. ";
0337         return vector<Trajectory>();
0338       }
0339 
0340       TrajectoryStateOnSurface smooTsos = combiner((*itm).updatedState(), predTsos);
0341 
0342       if (!smooTsos.isValid()) {
0343         LogTrace(category_) << "Error: smooth: combining smooth TSOS failed. ";
0344         return vector<Trajectory>();
0345       }
0346 
0347       myTraj.push(TrajectoryMeasurement((*itm).forwardPredictedState(),
0348                                         predTsos,
0349                                         smooTsos,
0350                                         (*itm).recHit(),
0351                                         theEstimator->estimate(combTsos, (*(*itm).recHit())).second  //,
0352                                         /*(*itm).layer()*/),
0353                   (*itm).estimate());
0354     } else {
0355       currTsos = predTsos;
0356       TrajectoryStateOnSurface combTsos = combiner(predTsos, (*itm).forwardPredictedState());
0357 
0358       if (!combTsos.isValid()) {
0359         LogTrace(category_) << "Error: smooth: combining TSOS failed. ";
0360         return vector<Trajectory>();
0361       }
0362 
0363       myTraj.push(TrajectoryMeasurement((*itm).forwardPredictedState(),
0364                                         predTsos,
0365                                         combTsos,
0366                                         (*itm).recHit()  //,
0367                                         /*(*itm).layer()*/));
0368     }
0369   }
0370 
0371   // last smoothed TrajectoryMeasurement is last filtered
0372   predTsos = propagatorOpposite()->propagate(currTsos, avtm.front().recHit()->det()->surface());
0373 
0374   if (!predTsos.isValid()) {
0375     LogTrace(category_) << "Error: last predict TSOS failed, use original one. ";
0376     //    return vector<Trajectory>();
0377     myTraj.push(TrajectoryMeasurement(avtm.front().forwardPredictedState(), avtm.front().recHit()));
0378   } else {
0379     if (avtm.front().recHit()->isValid()) {
0380       //update
0381       currTsos = theUpdator->update(predTsos, (*avtm.front().recHit()));
0382       if (currTsos.isValid())
0383         myTraj.push(TrajectoryMeasurement(avtm.front().forwardPredictedState(),
0384                                           predTsos,
0385                                           currTsos,
0386                                           avtm.front().recHit(),
0387                                           theEstimator->estimate(predTsos, (*avtm.front().recHit())).second  //,
0388                                           /*avtm.front().layer()*/),
0389                     avtm.front().estimate());
0390     }
0391   }
0392   LogTrace(category_) << "myTraj foundHits. " << myTraj.foundHits();
0393 
0394   if (myTraj.foundHits() >= 3)
0395     return vector<Trajectory>(1, myTraj);
0396   else {
0397     LogTrace(category_) << "Error: smooth: No enough hits in trajctory. ";
0398     return vector<Trajectory>();
0399   }
0400 }
0401 
0402 TrajectoryStateOnSurface CosmicMuonSmoother::initialState(const Trajectory& t) const {
0403   if (t.empty())
0404     return TrajectoryStateOnSurface();
0405 
0406   if (!t.firstMeasurement().updatedState().isValid() || !t.lastMeasurement().updatedState().isValid())
0407     return TrajectoryStateOnSurface();
0408 
0409   TrajectoryStateOnSurface result;
0410 
0411   bool beamhaloFlag = (t.firstMeasurement().updatedState().globalMomentum().eta() > 4.0 ||
0412                        t.lastMeasurement().updatedState().globalMomentum().eta() > 4.0);
0413 
0414   if (!beamhaloFlag) {  //initialState is the top one
0415     if (t.firstMeasurement().updatedState().globalPosition().y() >
0416         t.lastMeasurement().updatedState().globalPosition().y()) {
0417       result = t.firstMeasurement().updatedState();
0418     } else {
0419       result = t.lastMeasurement().updatedState();
0420     }
0421     if (result.globalMomentum().y() > 1.0)  //top tsos should pointing down
0422       theUtilities->reverseDirection(result, &*theService->magneticField());
0423   } else {
0424     if (t.firstMeasurement().updatedState().globalPosition().z() *
0425             t.lastMeasurement().updatedState().globalPosition().z() >
0426         0) {  //same side
0427       if (fabs(t.firstMeasurement().updatedState().globalPosition().z()) >
0428           fabs(t.lastMeasurement().updatedState().globalPosition().z())) {
0429         result = t.firstMeasurement().updatedState();
0430       } else {
0431         result = t.lastMeasurement().updatedState();
0432       }
0433     } else {  //different sides
0434 
0435       if (fabs(t.firstMeasurement().updatedState().globalPosition().eta()) >
0436           fabs(t.lastMeasurement().updatedState().globalPosition().eta())) {
0437         result = t.firstMeasurement().updatedState();
0438       } else {
0439         result = t.lastMeasurement().updatedState();
0440       }
0441     }
0442   }
0443 
0444   return result;
0445 }
0446 
0447 void CosmicMuonSmoother::sortHitsAlongMom(ConstRecHitContainer& hits, const TrajectoryStateOnSurface& tsos) const {
0448   if (hits.size() < 2)
0449     return;
0450   float dis1 = (hits.front()->globalPosition() - tsos.globalPosition()).mag();
0451   float dis2 = (hits.back()->globalPosition() - tsos.globalPosition()).mag();
0452 
0453   if (dis1 > dis2)
0454     std::reverse(hits.begin(), hits.end());
0455 
0456   return;
0457 }
[ 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!