00001
00002
00003
00004
00005
00015 #include "FWCore/Framework/interface/EDProducer.h"
00016 #include "FWCore/Framework/interface/Event.h"
00017 #include "FWCore/Framework/interface/EventSetup.h"
00018 #include "FWCore/Framework/interface/ESHandle.h"
00019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00020 #include "FWCore/Utilities/interface/InputTag.h"
00021
00022 #include "DataFormats/Math/interface/deltaR.h"
00023
00024 #include "DataFormats/MuonReco/interface/Muon.h"
00025 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
00026 #include "MagneticField/Engine/interface/MagneticField.h"
00027 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
00028
00029 #include "TrackingTools/TrackRefitter/interface/TrackTransformer.h"
00030 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00031 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00032 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00033 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00034 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
00035 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
00036 #include "TrackingTools/MeasurementDet/interface/MeasurementDet.h"
00037 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
00038 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
00039 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
00040 #include "TrackingTools/PatternTools/interface/TrajectoryStateUpdator.h"
00041 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimatorBase.h"
00042 #include "CommonTools/Utils/interface/StringCutObjectSelector.h"
00043 #include "TrackingTools/PatternTools/interface/TrajectoryMeasurement.h"
00044 #include "TrackingTools/PatternTools/interface/TrajMeasLessEstim.h"
00045
00046 class OutsideInMuonSeeder : public edm::EDProducer {
00047 public:
00048 explicit OutsideInMuonSeeder(const edm::ParameterSet & iConfig);
00049 virtual ~OutsideInMuonSeeder() { }
00050
00051 virtual void produce(edm::Event & iEvent, const edm::EventSetup & iSetup);
00052
00053 private:
00055 edm::InputTag src_;
00056
00058 StringCutObjectSelector<reco::Muon> selector_;
00059
00061 int layersToTry_;
00062
00064 int hitsToTry_;
00065
00067 bool fromVertex_;
00068
00070 double errorRescaling_;
00071
00072 std::string trackerPropagatorName_;
00073 std::string muonPropagatorName_;
00074 std::string measurementTrackerName_;
00075 std::string estimatorName_;
00076 std::string updatorName_;
00077
00078 double minEtaForTEC_, maxEtaForTOB_;
00079
00080 edm::ESHandle<MagneticField> magfield_;
00081 edm::ESHandle<Propagator> muonPropagator_;
00082 edm::ESHandle<Propagator> trackerPropagator_;
00083 edm::ESHandle<GlobalTrackingGeometry> geometry_;
00084 edm::ESHandle<MeasurementTracker> measurementTracker_;
00085 edm::ESHandle<Chi2MeasurementEstimatorBase> estimator_;
00086 edm::ESHandle<TrajectoryStateUpdator> updator_;
00087
00089 bool debug_;
00090
00092 Plane::PlanePointer dummyPlane_;
00093
00094 int doLayer(const GeometricSearchDet &layer, const TrajectoryStateOnSurface &state, std::vector<TrajectorySeed> &out) const ;
00095 void doDebug(const reco::Track &tk) const;
00096
00097 };
00098
00099 OutsideInMuonSeeder::OutsideInMuonSeeder(const edm::ParameterSet & iConfig) :
00100 src_(iConfig.getParameter<edm::InputTag>("src")),
00101 selector_(iConfig.existsAs<std::string>("cut") ? iConfig.getParameter<std::string>("cut") : "", true),
00102 layersToTry_(iConfig.getParameter<int32_t>("layersToTry")),
00103 hitsToTry_(iConfig.getParameter<int32_t>("hitsToTry")),
00104 fromVertex_(iConfig.getParameter<bool>("fromVertex")),
00105 errorRescaling_(iConfig.getParameter<double>("errorRescaleFactor")),
00106 trackerPropagatorName_(iConfig.getParameter<std::string>("trackerPropagator")),
00107 muonPropagatorName_(iConfig.getParameter<std::string>("muonPropagator")),
00108 estimatorName_(iConfig.getParameter<std::string>("hitCollector")),
00109 minEtaForTEC_(iConfig.getParameter<double>("minEtaForTEC")),
00110 maxEtaForTOB_(iConfig.getParameter<double>("maxEtaForTOB")),
00111 debug_(iConfig.getUntrackedParameter<bool>("debug",false)),
00112 dummyPlane_(Plane::build(Plane::PositionType(), Plane::RotationType()))
00113 {
00114 produces<std::vector<TrajectorySeed> >();
00115 measurementTrackerName_ = "";
00116 updatorName_ = "KFUpdator";
00117 }
00118
00119 void
00120 OutsideInMuonSeeder::produce(edm::Event & iEvent, const edm::EventSetup & iSetup) {
00121 using namespace edm;
00122 using namespace std;
00123 iSetup.get<IdealMagneticFieldRecord>().get(magfield_);
00124 iSetup.get<TrackingComponentsRecord>().get(trackerPropagatorName_, trackerPropagator_);
00125 iSetup.get<TrackingComponentsRecord>().get(muonPropagatorName_, muonPropagator_);
00126 iSetup.get<GlobalTrackingGeometryRecord>().get(geometry_);
00127 iSetup.get<CkfComponentsRecord>().get(measurementTracker_);
00128 iSetup.get<TrackingComponentsRecord>().get(estimatorName_,estimator_);
00129 iSetup.get<TrackingComponentsRecord>().get(updatorName_,updator_);
00130
00131 measurementTracker_->update(iEvent);
00132
00133 Handle<View<reco::Muon> > src;
00134 iEvent.getByLabel(src_, src);
00135
00136
00137 auto_ptr<vector<TrajectorySeed> > out(new vector<TrajectorySeed>());
00138
00139 for (View<reco::Muon>::const_iterator it = src->begin(), ed = src->end(); it != ed; ++it) {
00140 const reco::Muon &mu = *it;
00141 if (mu.outerTrack().isNull() || !selector_(mu)) continue;
00142 if (debug_ && mu.innerTrack().isNonnull()) doDebug(*mu.innerTrack());
00143
00144 muonPropagator_->setPropagationDirection(fromVertex_ ? alongMomentum : oppositeToMomentum);
00145 trackerPropagator_->setPropagationDirection(alongMomentum);
00146
00147 int sizeBefore = out->size();
00148 if (debug_) std::cout << "\n\n\nSeeding for muon of pt " << mu.pt() << ", eta " << mu.eta() << ", phi " << mu.phi() << std::endl;
00149 const reco::Track &tk = *mu.outerTrack();
00150
00151 TrajectoryStateOnSurface state;
00152 if (fromVertex_) {
00153 FreeTrajectoryState fstate = trajectoryStateTransform::initialFreeState(tk, magfield_.product());
00154 dummyPlane_->move(fstate.position() - dummyPlane_->position());
00155 state = TrajectoryStateOnSurface(fstate, *dummyPlane_);
00156 } else {
00157 state = trajectoryStateTransform::innerStateOnSurface(tk, *geometry_, magfield_.product());
00158 }
00159 if (std::abs(tk.eta()) < maxEtaForTOB_) {
00160 std::vector< BarrelDetLayer * > const & tob = measurementTracker_->geometricSearchTracker()->tobLayers();
00161 int iLayer = 6, found = 0;
00162 for (std::vector<BarrelDetLayer *>::const_reverse_iterator it = tob.rbegin(), ed = tob.rend(); it != ed; ++it, --iLayer) {
00163 if (debug_) std::cout << "\n ==== Trying TOB " << iLayer << " ====" << std::endl;
00164 if (doLayer(**it, state, *out)) {
00165 if (++found == layersToTry_) break;
00166 }
00167 }
00168 }
00169 if (tk.eta() > minEtaForTEC_) {
00170 int iLayer = 9, found = 0;
00171 std::vector< ForwardDetLayer * > const & tec = measurementTracker_->geometricSearchTracker()->posTecLayers();
00172 for (std::vector<ForwardDetLayer *>::const_reverse_iterator it = tec.rbegin(), ed = tec.rend(); it != ed; ++it, --iLayer) {
00173 if (debug_) std::cout << "\n ==== Trying TEC " << +iLayer << " ====" << std::endl;
00174 if (doLayer(**it, state, *out)) {
00175 if (++found == layersToTry_) break;
00176 }
00177 }
00178 }
00179 if (tk.eta() < -minEtaForTEC_) {
00180 int iLayer = 9, found = 0;
00181 std::vector< ForwardDetLayer * > const & tec = measurementTracker_->geometricSearchTracker()->negTecLayers();
00182 for (std::vector<ForwardDetLayer *>::const_reverse_iterator it = tec.rbegin(), ed = tec.rend(); it != ed; ++it, --iLayer) {
00183 if (debug_) std::cout << "\n ==== Trying TEC " << -iLayer << " ====" << std::endl;
00184 if (doLayer(**it, state, *out)) {
00185 if (++found == layersToTry_) break;
00186 }
00187 }
00188 }
00189 if (debug_) std::cout << "Outcome of seeding for muon of pt " << mu.pt() << ", eta " << mu.eta() << ", phi " << mu.phi() << ": found " << (out->size() - sizeBefore) << " seeds."<< std::endl;
00190
00191 }
00192
00193 iEvent.put(out);
00194 }
00195
00196 int
00197 OutsideInMuonSeeder::doLayer(const GeometricSearchDet &layer, const TrajectoryStateOnSurface &state, std::vector<TrajectorySeed> &out) const {
00198 TrajectoryStateOnSurface onLayer(state);
00199 onLayer.rescaleError(errorRescaling_);
00200 std::vector< GeometricSearchDet::DetWithState > dets;
00201 layer.compatibleDetsV(onLayer, *muonPropagator_, *estimator_, dets);
00202
00203 if (debug_) {
00204 std::cout << "Query on layer around x = " << onLayer.globalPosition() <<
00205 " with local pos error " << sqrt(onLayer.localError().positionError().xx()) << " , " << sqrt(onLayer.localError().positionError().yy()) << " , " <<
00206 " returned " << dets.size() << " compatible detectors" << std::endl;
00207 }
00208
00209 std::vector<TrajectoryMeasurement> meas;
00210 for (std::vector<GeometricSearchDet::DetWithState>::const_iterator it = dets.begin(), ed = dets.end(); it != ed; ++it) {
00211 const MeasurementDet *det = measurementTracker_->idToDet(it->first->geographicalId());
00212 if (det == 0) { std::cerr << "BOGUS detid " << it->first->geographicalId().rawId() << std::endl; continue; }
00213 if (!it->second.isValid()) continue;
00214 std::vector < TrajectoryMeasurement > mymeas = det->fastMeasurements(it->second, state, *trackerPropagator_, *estimator_);
00215 if (debug_) std::cout << "Query on detector " << it->first->geographicalId().rawId() << " returned " << mymeas.size() << " measurements." << std::endl;
00216 for (std::vector<TrajectoryMeasurement>::const_iterator it2 = mymeas.begin(), ed2 = mymeas.end(); it2 != ed2; ++it2) {
00217 if (it2->recHit()->isValid()) meas.push_back(*it2);
00218 }
00219 }
00220 int found = 0;
00221 std::sort(meas.begin(), meas.end(), TrajMeasLessEstim());
00222 for (std::vector<TrajectoryMeasurement>::const_iterator it2 = meas.begin(), ed2 = meas.end(); it2 != ed2; ++it2) {
00223 if (debug_) {
00224 std::cout << " inspecting Hit with chi2 = " << it2->estimate() << std::endl;
00225 std::cout << " track state " << it2->forwardPredictedState().globalPosition() << std::endl;
00226 std::cout << " rechit position " << it2->recHit()->globalPosition() << std::endl;
00227 }
00228 TrajectoryStateOnSurface updated = updator_->update(it2->forwardPredictedState(), *it2->recHit());
00229 if (updated.isValid()) {
00230 if (debug_) std::cout << " --> updated state: x = " << updated.globalPosition() << ", p = " << updated.globalMomentum() << std::endl;
00231 edm::OwnVector<TrackingRecHit> seedHits;
00232 seedHits.push_back(*it2->recHit()->hit());
00233 PTrajectoryStateOnDet const & pstate = trajectoryStateTransform::persistentState(updated, it2->recHit()->geographicalId().rawId());
00234 TrajectorySeed seed(pstate, std::move(seedHits), oppositeToMomentum);
00235 out.push_back(seed);
00236 found++; if (found == hitsToTry_) break;
00237 }
00238 }
00239 return found;
00240 }
00241
00242 void
00243 OutsideInMuonSeeder::doDebug(const reco::Track &tk) const {
00244 TrajectoryStateOnSurface tsos = trajectoryStateTransform::innerStateOnSurface(tk, *geometry_, &*magfield_);
00245 muonPropagator_->setPropagationDirection(alongMomentum);
00246 for (unsigned int i = 0; i < tk.recHitsSize(); ++i) {
00247 const TrackingRecHit *hit = &*tk.recHit(i);
00248 const GeomDet *det = geometry_->idToDet(hit->geographicalId());
00249 if (det == 0) continue;
00250 if (i != 0) tsos = muonPropagator_->propagate(tsos, det->surface());
00251 if (!tsos.isValid()) continue;
00252 std::cout << " state " << i << " at x = " << tsos.globalPosition() << ", p = " << tsos.globalMomentum() << std::endl;
00253 if (hit->isValid()) {
00254 std::cout << " valid rechit on detid " << hit->geographicalId().rawId() << std::endl;
00255 } else {
00256 std::cout << " invalid rechit on detid " << hit->geographicalId().rawId() << std::endl;
00257 }
00258 }
00259 }
00260
00261 #include "FWCore/Framework/interface/MakerMacros.h"
00262 DEFINE_FWK_MODULE(OutsideInMuonSeeder);
1.7.3