L3TrackFinder/src/MuonCkfTrajectoryBuilder.cc

0001 #include "RecoMuon/L3TrackFinder/interface/MuonCkfTrajectoryBuilder.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0005 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0006 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
0007 #include "TrackingTools/PatternTools/interface/TrajMeasLessEstim.h"
0008 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
0009 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0010 #include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h"
0011 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimatorBase.h"
0012 #include "TrackingTools/TrajectoryFiltering/interface/TrajectoryFilter.h"
0013 #include "TrackingTools/PatternTools/interface/TransverseImpactPointExtrapolator.h"
0014 #include "TrackingTools/DetLayers/interface/NavigationSchool.h"
0015 #include "RecoMuon/L3TrackFinder/src/EtaPhiEstimator.h"
0016 #include <sstream>
0017
0018 MuonCkfTrajectoryBuilder::MuonCkfTrajectoryBuilder(const edm::ParameterSet& conf, edm::ConsumesCollector& iC)
0019     : CkfTrajectoryBuilder(conf, iC),
0020       theDeltaEta(conf.getParameter<double>("deltaEta")),
0021       theDeltaPhi(conf.getParameter<double>("deltaPhi")),
0022       theProximityPropagatorName(conf.getParameter<std::string>("propagatorProximity")),
0023       theProximityPropagator(nullptr),
0024       thePropagatorToken(iC.esConsumes(edm::ESInputTag("", theProximityPropagatorName))),
0025       theEtaPhiEstimator(nullptr) {
0026   //and something specific to me ?
0027   theUseSeedLayer = conf.getParameter<bool>("useSeedLayer");
0028   theRescaleErrorIfFail = conf.getParameter<double>("rescaleErrorIfFail");
0029 }
0030
0031 MuonCkfTrajectoryBuilder::~MuonCkfTrajectoryBuilder() {}
0032
0033 void MuonCkfTrajectoryBuilder::fillPSetDescription(edm::ParameterSetDescription& iDesc) {
0034   CkfTrajectoryBuilder::fillPSetDescription(iDesc);
0035   iDesc.add<double>("deltaEta", .1);
0036   iDesc.add<double>("deltaPhi", .1);
0037   iDesc.add<std::string>("propagatorProximity", "SteppingHelixPropagatorAny");
0038   iDesc.add<bool>("useSeedLayer", false);
0039   iDesc.add<double>("rescaleErrorIfFail", 1.);
0040 }
0041
0042 void MuonCkfTrajectoryBuilder::setEvent_(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0043   CkfTrajectoryBuilder::setEvent_(iEvent, iSetup);
0044
0045   theProximityPropagator = &iSetup.getData(thePropagatorToken);
0046
0047   // theEstimator is set for this event in the base class
0048   if (theEstimatorWatcher.check(iSetup) && theDeltaEta > 0 && theDeltaPhi > 0)
0049     theEtaPhiEstimator = std::make_unique<EtaPhiEstimator>(theDeltaEta, theDeltaPhi, theEstimator);
0050 }
0051
0052 /*
0053 std::string dumpMeasurement(const TrajectoryMeasurement & tm)
0054 {
0055   std::stringstream buffer;
0056   buffer<<"layer pointer: "<<tm.layer()<<"\n"
0057         <<"estimate: "<<tm.estimate()<<"\n"
0058         <<"forward state: \n"
0059         <<"x: "<<tm.forwardPredictedState().globalPosition()<<"\n"
0060         <<"p: "<<tm.forwardPredictedState().globalMomentum()<<"\n"
0061         <<"geomdet pointer from rechit: "<<tm.recHit()->det()<<"\n"
0062         <<"detId: "<<tm.recHit()->geographicalId().rawId();
0063   if (tm.recHit()->isValid()){
0064     buffer<<"\n hit global x: "<<tm.recHit()->globalPosition()
0065       <<"\n hit global error: "<<tm.recHit()->globalPositionError().matrix()
0066       <<"\n hit local x:"<<tm.recHit()->localPosition()
0067       <<"\n hit local error"<<tm.recHit()->localPositionError();
0068   }
0069   return buffer.str();
0070 }
0071 std::string dumpMeasurements(const std::vector<TrajectoryMeasurement> & v)
0072 {
0073   std::stringstream buffer;
0074   buffer<<v.size()<<" total measurements\n";
0075   for (std::vector<TrajectoryMeasurement>::const_iterator it = v.begin(); it!=v.end();++it){
0076     buffer<<dumpMeasurement(*it);
0077     buffer<<"\n";}
0078   return buffer.str();
0079 }
0080 */
0081
0082 void MuonCkfTrajectoryBuilder::collectMeasurement(const DetLayer* layer,
0083                                                   const std::vector<const DetLayer*>& nl,
0084                                                   const TrajectoryStateOnSurface& currentState,
0085                                                   std::vector<TM>& result,
0086                                                   int& invalidHits,
0087                                                   const Propagator* prop) const {
0088   for (std::vector<const DetLayer*>::const_iterator il = nl.begin(); il != nl.end(); il++) {
0089     TSOS stateToUse = currentState;
0090
0091     if (layer == (*il)) {
0092       LogDebug("CkfPattern") << " self propagating in findCompatibleMeasurements.\n from: \n" << stateToUse;
0093       //self navigation case
0094       // go to a middle point first
0095       TransverseImpactPointExtrapolator middle;
0096       GlobalPoint center(0, 0, 0);
0097       stateToUse = middle.extrapolate(stateToUse, center, *prop);
0098
0099       if (!stateToUse.isValid())
0100         continue;
0101       LogDebug("CkfPattern") << "to: " << stateToUse;
0102     }
0103
0104     LayerMeasurements layerMeasurements(theMeasurementTracker->measurementTracker(), *theMeasurementTracker);
0105     std::vector<TM> tmp = layerMeasurements.measurements((**il), stateToUse, *prop, *theEstimator);
0106
0107     if (tmp.size() == 1 && theEtaPhiEstimator) {
0108       LogDebug("CkfPattern") << "only an invalid hit is found. trying differently";
0109       tmp = layerMeasurements.measurements((**il), stateToUse, *prop, *theEtaPhiEstimator);
0110     }
0111     LogDebug("CkfPattern") << tmp.size() << " measurements returned by LayerMeasurements";
0112
0113     if (!tmp.empty()) {
0114       // FIXME durty-durty-durty cleaning: never do that please !
0115       /*      for (vector<TM>::iterator it = tmp.begin(); it!=tmp.end(); ++it)
0116               {if (it->recHit()->det()==0) it=tmp.erase(it)--;}*/
0117
0118       if (result.empty())
0119         result = tmp;
0120       else {
0121         // keep one dummy TM at the end, skip the others
0122         result.insert(result.end() - invalidHits, tmp.begin(), tmp.end());
0123       }
0124       invalidHits++;
0125     }
0126   }
0127
0128   LogDebug("CkfPattern") << "starting from:\n"
0129                          << "x: " << currentState.globalPosition() << "\n"
0130                          << "p: " << currentState.globalMomentum() << "\n"
0131                          << PrintoutHelper::dumpMeasurements(result);
0132 }
0133
0134 void MuonCkfTrajectoryBuilder::findCompatibleMeasurements(const TrajectorySeed& seed,
0135                                                           const TempTrajectory& traj,
0136                                                           std::vector<TrajectoryMeasurement>& result) const {
0137   int invalidHits = 0;
0138
0139   std::vector<const DetLayer*> nl;
0140
0141   if (traj.empty()) {
0142     LogDebug("CkfPattern") << "using JR patch for no measurement case";
0143     //what if there are no measurement on the Trajectory
0144
0145     //set the currentState to be the one from the trajectory seed starting point
0146     PTrajectoryStateOnDet ptod = seed.startingState();
0147     DetId id(ptod.detId());
0148     const GeomDet* g = theMeasurementTracker->geomTracker()->idToDet(id);
0149     const Surface* surface = &g->surface();
0150
0151     TrajectoryStateOnSurface currentState(
0152         trajectoryStateTransform::transientState(ptod, surface, forwardPropagator(seed)->magneticField()));
0153
0154     //set the next layers to be that one the state is on
0155     const DetLayer* l = theMeasurementTracker->geometricSearchTracker()->detLayer(id);
0156
0157     if (theUseSeedLayer) {
0158       {
0159         //get the measurements on the layer first
0160         LogDebug("CkfPattern") << "using the layer of the seed first.";
0161         nl.push_back(l);
0162         collectMeasurement(l, nl, currentState, result, invalidHits, theProximityPropagator);
0163       }
0164
0165       //if fails: try to rescale locally the state to find measurements
0166       if ((unsigned int)invalidHits == result.size() && theRescaleErrorIfFail != 1.0 && !result.empty()) {
0167         result.clear();
0168         LogDebug("CkfPattern") << "using a rescale by " << theRescaleErrorIfFail << " to find measurements.";
0169         TrajectoryStateOnSurface rescaledCurrentState = currentState;
0170         rescaledCurrentState.rescaleError(theRescaleErrorIfFail);
0171         invalidHits = 0;
0172         collectMeasurement(l, nl, rescaledCurrentState, result, invalidHits, theProximityPropagator);
0173       }
0174     }
0175
0176     //if fails: go to next layers
0177     if (result.empty() || (unsigned int)invalidHits == result.size()) {
0178       result.clear();
0179       LogDebug("CkfPattern") << "Need to go to next layer to get measurements";
0180       //the following will "JUMP" the first layer measurements
0181       nl = theNavigationSchool->nextLayers(*l, *currentState.freeState(), traj.direction());
0182       if (nl.empty()) {
0183         LogDebug("CkfPattern") << " there was no next layer with wellInside. Use the next with no check.";
0184         //means you did not get any compatible layer on the next 1/2 tracker layer.
0185         // use the next layers with no checking
0186         nl = theNavigationSchool->nextLayers(*l, ((traj.direction() == alongMomentum) ? insideOut : outsideIn));
0187       }
0188       invalidHits = 0;
0189       collectMeasurement(l, nl, currentState, result, invalidHits, forwardPropagator(seed));
0190     }
0191
0192     //if fails: this is on the next layers already, try rescaling locally the state
0193     if (!result.empty() && (unsigned int)invalidHits == result.size() && theRescaleErrorIfFail != 1.0) {
0194       result.clear();
0195       LogDebug("CkfPattern") << "using a rescale by " << theRescaleErrorIfFail
0196                              << " to find measurements on next layers.";
0197       TrajectoryStateOnSurface rescaledCurrentState = currentState;
0198       rescaledCurrentState.rescaleError(theRescaleErrorIfFail);
0199       invalidHits = 0;
0200       collectMeasurement(l, nl, rescaledCurrentState, result, invalidHits, forwardPropagator(seed));
0201     }
0202
0203   } else  //regular case
0204   {
0205     TSOS currentState(traj.lastMeasurement().updatedState());
0206
0207     nl = theNavigationSchool->nextLayers(*traj.lastLayer(), *currentState.freeState(), traj.direction());
0208     if (nl.empty()) {
0209       LogDebug("CkfPattern") << " no next layers... going " << traj.direction() << "\n from: \n"
0210                              << currentState
0211                              << "\n from detId: " << traj.lastMeasurement().recHit()->geographicalId().rawId();
0212       return;
0213     }
0214
0215     collectMeasurement(traj.lastLayer(), nl, currentState, result, invalidHits, forwardPropagator(seed));
0216   }
0217
0218   // sort the final result, keep dummy measurements at the end
0219   if (result.size() > 1) {
0220     sort(result.begin(), result.end() - invalidHits, TrajMeasLessEstim());
0221   }
0222
0223 #ifdef DEBUG_INVALID
0224   bool afterInvalid = false;
0225   for (std::vector<TM>::const_iterator i = result.begin(); i != result.end(); i++) {
0226     if (!i->recHit().isValid())
0227       afterInvalid = true;
0228     if (afterInvalid && i->recHit().isValid()) {
0229       edm::LogError("CkfPattern") << "CkfTrajectoryBuilder error: valid hit after invalid!";
0230     }
0231   }
0232 #endif
0233
0234   //analyseMeasurements( result, traj);
0235 }