Skip to content

Commit

Permalink
Merge pull request #3037 from rovere/threadSafePropagators
Browse files Browse the repository at this point in the history
Multithreading -- Thread safe propagators
  • Loading branch information
ktf committed Mar 31, 2014
2 parents 4d9e8fe + 96becc8 commit b2b753a
Show file tree
Hide file tree
Showing 14 changed files with 359 additions and 329 deletions.
89 changes: 56 additions & 33 deletions RecoTracker/SpecialSeedGenerators/src/OutsideInMuonSeeder.cc
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,8 @@

/**
\class pat::OutsideInMuonSeeder MuonReSeeder.h "MuonAnalysis/MuonAssociators/interface/MuonReSeeder.h"
\brief Matcher of reconstructed objects to other reconstructed objects using the tracks inside them
\brief Matcher of reconstructed objects to other reconstructed objects using the tracks inside them
\author Giovanni Petrucciani
*/

Expand Down Expand Up @@ -68,11 +68,11 @@ class OutsideInMuonSeeder : public edm::EDProducer {
double errorRescaling_;

std::string trackerPropagatorName_;
std::string muonPropagatorName_;
std::string muonPropagatorName_;
edm::EDGetTokenT<MeasurementTrackerEvent> measurementTrackerTag_;
std::string measurementTrackerName_;
std::string estimatorName_;
std::string updatorName_;
std::string measurementTrackerName_;
std::string estimatorName_;
std::string updatorName_;

double minEtaForTEC_, maxEtaForTOB_;

Expand All @@ -85,12 +85,17 @@ class OutsideInMuonSeeder : public edm::EDProducer {

/// Dump deug information
bool debug_;

/// Surface used to make a TSOS at the PCA to the beamline
Plane::PlanePointer dummyPlane_;

int doLayer(const GeometricSearchDet &layer, const TrajectoryStateOnSurface &state, std::vector<TrajectorySeed> &out, const MeasurementTrackerEvent &mte) const ;
void doDebug(const reco::Track &tk) const;
int doLayer(const GeometricSearchDet &layer,
const TrajectoryStateOnSurface &state,
std::vector<TrajectorySeed> &out,
const Propagator &muon_propagator,
const Propagator &tracker_propagator,
const MeasurementTrackerEvent &mte) const ;
void doDebug(const reco::Track &tk) const;

};

Expand All @@ -110,11 +115,11 @@ OutsideInMuonSeeder::OutsideInMuonSeeder(const edm::ParameterSet & iConfig) :
debug_(iConfig.getUntrackedParameter<bool>("debug",false)),
dummyPlane_(Plane::build(Plane::PositionType(), Plane::RotationType()))
{
produces<std::vector<TrajectorySeed> >();
produces<std::vector<TrajectorySeed> >();
updatorName_ = "KFUpdator";
}

void
void
OutsideInMuonSeeder::produce(edm::Event & iEvent, const edm::EventSetup & iSetup) {
using namespace edm;
using namespace std;
Expand All @@ -123,30 +128,34 @@ OutsideInMuonSeeder::produce(edm::Event & iEvent, const edm::EventSetup & iSetup
iSetup.get<TrackingComponentsRecord>().get(trackerPropagatorName_, trackerPropagator_);
iSetup.get<TrackingComponentsRecord>().get(muonPropagatorName_, muonPropagator_);
iSetup.get<GlobalTrackingGeometryRecord>().get(geometry_);
iSetup.get<TrackingComponentsRecord>().get(estimatorName_,estimator_);
iSetup.get<TrackingComponentsRecord>().get(updatorName_,updator_);
iSetup.get<TrackingComponentsRecord>().get(estimatorName_,estimator_);
iSetup.get<TrackingComponentsRecord>().get(updatorName_,updator_);

Handle<MeasurementTrackerEvent> measurementTracker;
iEvent.getByToken(measurementTrackerTag_, measurementTracker);

Handle<View<reco::Muon> > src;
iEvent.getByToken(src_, src);


auto_ptr<vector<TrajectorySeed> > out(new vector<TrajectorySeed>());

for (View<reco::Muon>::const_iterator it = src->begin(), ed = src->end(); it != ed; ++it) {
const reco::Muon &mu = *it;
if (mu.outerTrack().isNull() || !selector_(mu)) continue;
if (debug_ && mu.innerTrack().isNonnull()) doDebug(*mu.innerTrack());

muonPropagator_->setPropagationDirection(fromVertex_ ? alongMomentum : oppositeToMomentum);
trackerPropagator_->setPropagationDirection(alongMomentum);
// Better clone here and not directly into doLayer to avoid
// useless clone/destroy operations to set, in the end, the
// very same direction every single time.
std::unique_ptr<Propagator> pmuon_cloned = SetPropagationDirection(*muonPropagator_,
fromVertex_ ? alongMomentum : oppositeToMomentum);
std::unique_ptr<Propagator> ptracker_cloned = SetPropagationDirection(*trackerPropagator_, alongMomentum);

int sizeBefore = out->size();
if (debug_) std::cout << "\n\n\nSeeding for muon of pt " << mu.pt() << ", eta " << mu.eta() << ", phi " << mu.phi() << std::endl;
const reco::Track &tk = *mu.outerTrack();

TrajectoryStateOnSurface state;
if (fromVertex_) {
FreeTrajectoryState fstate = trajectoryStateTransform::initialFreeState(tk, magfield_.product());
Expand All @@ -160,7 +169,10 @@ OutsideInMuonSeeder::produce(edm::Event & iEvent, const edm::EventSetup & iSetup
int iLayer = 6, found = 0;
for (std::vector<BarrelDetLayer *>::const_reverse_iterator it = tob.rbegin(), ed = tob.rend(); it != ed; ++it, --iLayer) {
if (debug_) std::cout << "\n ==== Trying TOB " << iLayer << " ====" << std::endl;
if (doLayer(**it, state, *out, *measurementTracker)) {
if (doLayer(**it, state, *out,
*(pmuon_cloned.get()),
*(ptracker_cloned.get()),
*measurementTracker)) {
if (++found == layersToTry_) break;
}
}
Expand All @@ -170,7 +182,10 @@ OutsideInMuonSeeder::produce(edm::Event & iEvent, const edm::EventSetup & iSetup
std::vector< ForwardDetLayer * > const & tec = measurementTracker->geometricSearchTracker()->posTecLayers();
for (std::vector<ForwardDetLayer *>::const_reverse_iterator it = tec.rbegin(), ed = tec.rend(); it != ed; ++it, --iLayer) {
if (debug_) std::cout << "\n ==== Trying TEC " << +iLayer << " ====" << std::endl;
if (doLayer(**it, state, *out, *measurementTracker)) {
if (doLayer(**it, state, *out,
*(pmuon_cloned.get()),
*(ptracker_cloned.get()),
*measurementTracker)) {
if (++found == layersToTry_) break;
}
}
Expand All @@ -180,28 +195,36 @@ OutsideInMuonSeeder::produce(edm::Event & iEvent, const edm::EventSetup & iSetup
std::vector< ForwardDetLayer * > const & tec = measurementTracker->geometricSearchTracker()->negTecLayers();
for (std::vector<ForwardDetLayer *>::const_reverse_iterator it = tec.rbegin(), ed = tec.rend(); it != ed; ++it, --iLayer) {
if (debug_) std::cout << "\n ==== Trying TEC " << -iLayer << " ====" << std::endl;
if (doLayer(**it, state, *out, *measurementTracker)) {
if (doLayer(**it, state, *out,
*(pmuon_cloned.get()),
*(ptracker_cloned.get()),
*measurementTracker)) {
if (++found == layersToTry_) break;
}
}
}
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;

}

iEvent.put(out);
}

int
OutsideInMuonSeeder::doLayer(const GeometricSearchDet &layer, const TrajectoryStateOnSurface &state, std::vector<TrajectorySeed> &out, const MeasurementTrackerEvent &measurementTracker) const {
OutsideInMuonSeeder::doLayer(const GeometricSearchDet &layer,
const TrajectoryStateOnSurface &state,
std::vector<TrajectorySeed> &out,
const Propagator & muon_propagator,
const Propagator & tracker_propagator,
const MeasurementTrackerEvent &measurementTracker) const {
TrajectoryStateOnSurface onLayer(state);
onLayer.rescaleError(errorRescaling_);
std::vector< GeometricSearchDet::DetWithState > dets;
layer.compatibleDetsV(onLayer, *muonPropagator_, *estimator_, dets);
layer.compatibleDetsV(onLayer, muon_propagator, *estimator_, dets);

if (debug_) {
std::cout << "Query on layer around x = " << onLayer.globalPosition() <<
" with local pos error " << sqrt(onLayer.localError().positionError().xx()) << " , " << sqrt(onLayer.localError().positionError().yy()) << " , " <<
std::cout << "Query on layer around x = " << onLayer.globalPosition() <<
" with local pos error " << sqrt(onLayer.localError().positionError().xx()) << " , " << sqrt(onLayer.localError().positionError().yy()) << " , " <<
" returned " << dets.size() << " compatible detectors" << std::endl;
}

Expand All @@ -210,7 +233,7 @@ OutsideInMuonSeeder::doLayer(const GeometricSearchDet &layer, const TrajectorySt
MeasurementDetWithData det = measurementTracker.idToDet(it->first->geographicalId());
if (det.isNull()) { std::cerr << "BOGUS detid " << it->first->geographicalId().rawId() << std::endl; continue; }
if (!it->second.isValid()) continue;
std::vector < TrajectoryMeasurement > mymeas = det.fastMeasurements(it->second, state, *trackerPropagator_, *estimator_);
std::vector < TrajectoryMeasurement > mymeas = det.fastMeasurements(it->second, state, tracker_propagator, *estimator_);
if (debug_) std::cout << "Query on detector " << it->first->geographicalId().rawId() << " returned " << mymeas.size() << " measurements." << std::endl;
for (std::vector<TrajectoryMeasurement>::const_iterator it2 = mymeas.begin(), ed2 = mymeas.end(); it2 != ed2; ++it2) {
if (it2->recHit()->isValid()) meas.push_back(*it2);
Expand All @@ -221,16 +244,16 @@ OutsideInMuonSeeder::doLayer(const GeometricSearchDet &layer, const TrajectorySt
for (std::vector<TrajectoryMeasurement>::const_iterator it2 = meas.begin(), ed2 = meas.end(); it2 != ed2; ++it2) {
if (debug_) {
std::cout << " inspecting Hit with chi2 = " << it2->estimate() << std::endl;
std::cout << " track state " << it2->forwardPredictedState().globalPosition() << std::endl;
std::cout << " rechit position " << it2->recHit()->globalPosition() << std::endl;
std::cout << " track state " << it2->forwardPredictedState().globalPosition() << std::endl;
std::cout << " rechit position " << it2->recHit()->globalPosition() << std::endl;
}
TrajectoryStateOnSurface updated = updator_->update(it2->forwardPredictedState(), *it2->recHit());
if (updated.isValid()) {
if (debug_) std::cout << " --> updated state: x = " << updated.globalPosition() << ", p = " << updated.globalMomentum() << std::endl;
if (debug_) std::cout << " --> updated state: x = " << updated.globalPosition() << ", p = " << updated.globalMomentum() << std::endl;
edm::OwnVector<TrackingRecHit> seedHits;
seedHits.push_back(*it2->recHit()->hit());
PTrajectoryStateOnDet const & pstate = trajectoryStateTransform::persistentState(updated, it2->recHit()->geographicalId().rawId());
TrajectorySeed seed(pstate, std::move(seedHits), oppositeToMomentum);
TrajectorySeed seed(pstate, std::move(seedHits), oppositeToMomentum);
out.push_back(seed);
found++; if (found == hitsToTry_) break;
}
Expand All @@ -241,12 +264,12 @@ OutsideInMuonSeeder::doLayer(const GeometricSearchDet &layer, const TrajectorySt
void
OutsideInMuonSeeder::doDebug(const reco::Track &tk) const {
TrajectoryStateOnSurface tsos = trajectoryStateTransform::innerStateOnSurface(tk, *geometry_, &*magfield_);
muonPropagator_->setPropagationDirection(alongMomentum);
std::unique_ptr<Propagator> pmuon_cloned = SetPropagationDirection(*muonPropagator_, alongMomentum);
for (unsigned int i = 0; i < tk.recHitsSize(); ++i) {
const TrackingRecHit *hit = &*tk.recHit(i);
const GeomDet *det = geometry_->idToDet(hit->geographicalId());
const GeomDet *det = geometry_->idToDet(hit->geographicalId());
if (det == 0) continue;
if (i != 0) tsos = muonPropagator_->propagate(tsos, det->surface());
if (i != 0) tsos = pmuon_cloned->propagate(tsos, det->surface());
if (!tsos.isValid()) continue;
std::cout << " state " << i << " at x = " << tsos.globalPosition() << ", p = " << tsos.globalMomentum() << std::endl;
if (hit->isValid()) {
Expand Down
52 changes: 26 additions & 26 deletions TrackingTools/GeomPropagators/interface/BeamHaloPropagator.h
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@

/** \class BeamHaloPropagator
*
* A propagator which use different algorithm to propagate
* A propagator which use different algorithm to propagate
* within an endcap or to cross over to the other endcap
*
* \author Jean-Roch VLIMANT UCSB
Expand All @@ -25,7 +25,7 @@ class BeamHaloPropagator GCC11_FINAL : public Propagator {

public:

/* Constructor */
/* Constructor */
///Defines which propagator is used inside endcap and in barrel
BeamHaloPropagator(const Propagator* aEndCapTkProp, const Propagator* aCrossTkProp, const MagneticField* field,
PropagationDirection dir = alongMomentum);
Expand All @@ -38,7 +38,7 @@ class BeamHaloPropagator GCC11_FINAL : public Propagator {
///Copy constructor
BeamHaloPropagator( const BeamHaloPropagator& );

/** virtual destructor */
/** virtual destructor */
virtual ~BeamHaloPropagator() ;

///Virtual constructor (using copy c'tor)
Expand All @@ -47,20 +47,20 @@ class BeamHaloPropagator GCC11_FINAL : public Propagator {
}


void setPropagationDirection (PropagationDirection dir) const
void setPropagationDirection (PropagationDirection dir) override
{
Propagator::setPropagationDirection(dir);
getEndCapTkPropagator()->setPropagationDirection(dir);
getCrossTkPropagator()->setPropagationDirection(dir);
theEndCapTkProp->setPropagationDirection(dir);
theCrossTkProp->setPropagationDirection(dir);
}



/* Operations as propagator*/
TrajectoryStateOnSurface propagate(const FreeTrajectoryState& fts,
/* Operations as propagator*/
TrajectoryStateOnSurface propagate(const FreeTrajectoryState& fts,
const Surface& surface) const;

TrajectoryStateOnSurface propagate(const TrajectoryStateOnSurface& tsos,
TrajectoryStateOnSurface propagate(const TrajectoryStateOnSurface& tsos,
const Surface& surface) const {
return Propagator::propagate(tsos,surface);
}
Expand All @@ -73,42 +73,42 @@ class BeamHaloPropagator GCC11_FINAL : public Propagator {
return Propagator::propagate(tsos, plane);
}

TrajectoryStateOnSurface propagate(const FreeTrajectoryState& fts,
TrajectoryStateOnSurface propagate(const FreeTrajectoryState& fts,
const Cylinder& cylinder) const;

TrajectoryStateOnSurface propagate(const TrajectoryStateOnSurface& tsos,
TrajectoryStateOnSurface propagate(const TrajectoryStateOnSurface& tsos,
const Cylinder& cylinder) const {
return Propagator::propagate(tsos, cylinder);
}

std::pair<TrajectoryStateOnSurface,double>
propagateWithPath(const FreeTrajectoryState& fts,
std::pair<TrajectoryStateOnSurface,double>
propagateWithPath(const FreeTrajectoryState& fts,
const Surface& surface) const {
return Propagator::propagateWithPath(fts,surface);
}

std::pair<TrajectoryStateOnSurface,double>
propagateWithPath(const TrajectoryStateOnSurface& tsos,
std::pair<TrajectoryStateOnSurface,double>
propagateWithPath(const TrajectoryStateOnSurface& tsos,
const Surface& surface) const {
return Propagator::propagateWithPath(tsos,surface);
}

std::pair<TrajectoryStateOnSurface,double>
propagateWithPath(const FreeTrajectoryState& fts,
std::pair<TrajectoryStateOnSurface,double>
propagateWithPath(const FreeTrajectoryState& fts,
const Plane& plane) const;

std::pair<TrajectoryStateOnSurface,double>
propagateWithPath(const TrajectoryStateOnSurface& tsos,
std::pair<TrajectoryStateOnSurface,double>
propagateWithPath(const TrajectoryStateOnSurface& tsos,
const Plane& plane) const {
return Propagator::propagateWithPath(tsos, plane);
}

std::pair<TrajectoryStateOnSurface,double>
propagateWithPath(const FreeTrajectoryState& fts,
std::pair<TrajectoryStateOnSurface,double>
propagateWithPath(const FreeTrajectoryState& fts,
const Cylinder& cylinder) const;

std::pair<TrajectoryStateOnSurface,double>
propagateWithPath(const TrajectoryStateOnSurface& tsos,
std::pair<TrajectoryStateOnSurface,double>
propagateWithPath(const TrajectoryStateOnSurface& tsos,
const Cylinder& cylinder) const {
return Propagator::propagateWithPath(tsos, cylinder);
}
Expand All @@ -124,16 +124,16 @@ class BeamHaloPropagator GCC11_FINAL : public Propagator {
virtual const MagneticField* magneticField() const {return theField;}

private:
void directionCheck(PropagationDirection dir)const;
void directionCheck(PropagationDirection dir);

Propagator* theEndCapTkProp;
Propagator* theCrossTkProp;
const MagneticField* theField;

protected:

};

#endif
#endif


0 comments on commit b2b753a

Please sign in to comment.