Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2025-01-09 23:33:54

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