Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix in TSGForOI muons - 93X #19948

Merged
merged 3 commits into from Aug 2, 2017
Merged

Fix in TSGForOI muons - 93X #19948

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
e04a72c
fix TSGForOI
silviodonato Jul 28, 2017
5e45059
move TSGForOI to ::global
silviodonato Jul 31, 2017
644c0b8
improve readability
silviodonato Aug 1, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
148 changes: 89 additions & 59 deletions RecoMuon/TrackerSeedGenerator/plugins/TSGForOI.cc
View file
Open in desktop
Expand Up @@ -14,6 +14,7 @@ using namespace std;

TSGForOI::TSGForOI(const edm::ParameterSet & iConfig) :
src_(consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("src"))),
numOfMaxSeedsParam_(iConfig.getParameter<uint32_t>("maxSeeds")),
numOfLayersToTry_(iConfig.getParameter<int32_t>("layersToTry")),
numOfHitsToTry_(iConfig.getParameter<int32_t>("hitsToTry")),
fixedErrorRescalingForHits_(iConfig.getParameter<double>("fixedErrorRescaleFactorForHits")),
Expand All @@ -25,7 +26,6 @@ TSGForOI::TSGForOI(const edm::ParameterSet & iConfig) :
maxEtaForTOB_(iConfig.getParameter<double>("maxEtaForTOB")),
useHitLessSeeds_(iConfig.getParameter<bool>("UseHitLessSeeds")),
useStereoLayersInTEC_(iConfig.getParameter<bool>("UseStereoLayersInTEC")),
dummyPlane_(Plane::build(Plane::PositionType(), Plane::RotationType())),
updator_(new KFUpdator()),
measurementTrackerTag_(consumes<MeasurementTrackerEvent>(iConfig.getParameter<edm::InputTag>("MeasurementTrackerEvent"))),
pT1_(iConfig.getParameter<double>("pT1")),
Expand All @@ -38,43 +38,62 @@ TSGForOI::TSGForOI(const edm::ParameterSet & iConfig) :
SF3_(iConfig.getParameter<double>("SF3")),
SF4_(iConfig.getParameter<double>("SF4")),
SF5_(iConfig.getParameter<double>("SF5")),
tsosDiff_(iConfig.getParameter<double>("tsosDiff"))
tsosDiff_(iConfig.getParameter<double>("tsosDiff")),
theCategory(string("Muon|RecoMuon|TSGForOI"))
{
numOfMaxSeeds_=iConfig.getParameter<uint32_t>("maxSeeds");
produces<std::vector<TrajectorySeed> >();
numSeedsMade_=0;
theCategory = "Muon|RecoMuon|TSGForOI";
}


TSGForOI::~TSGForOI(){
}


void TSGForOI::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
iSetup.get<IdealMagneticFieldRecord>().get(magfield_);
iSetup.get<TrackingComponentsRecord>().get("PropagatorWithMaterial", propagatorOpposite_);
iSetup.get<TrackingComponentsRecord>().get("PropagatorWithMaterial", propagatorAlong_);
iSetup.get<GlobalTrackingGeometryRecord>().get(geometry_);
iSetup.get<TrackingComponentsRecord>().get(estimatorName_,estimator_);
iEvent.getByToken(measurementTrackerTag_, measurementTracker_);
edm::Handle<reco::TrackCollection> l2TrackCol;
void TSGForOI::produce(edm::StreamID sid, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
/// Init variables
unsigned int numOfMaxSeeds = numOfMaxSeedsParam_;
unsigned int numSeedsMade=0;
bool analysedL2 = false;
bool foundHitlessSeed = false;
unsigned int layerCount = 0;


/// Surface used to make a TSOS at the PCA to the beamline
Plane::PlanePointer dummyPlane = Plane::build(Plane::PositionType(), Plane::RotationType());

/// Read ESHandles
edm::Handle<MeasurementTrackerEvent> measurementTrackerH;
edm::ESHandle<Chi2MeasurementEstimatorBase> estimatorH;
edm::ESHandle<MagneticField> magfieldH;
edm::ESHandle<Propagator> propagatorAlongH;
edm::ESHandle<Propagator> propagatorOppositeH;
edm::ESHandle<GlobalTrackingGeometry> geometryH;

iSetup.get<IdealMagneticFieldRecord>().get(magfieldH);
iSetup.get<TrackingComponentsRecord>().get("PropagatorWithMaterial", propagatorOppositeH);
iSetup.get<TrackingComponentsRecord>().get("PropagatorWithMaterial", propagatorAlongH);
iSetup.get<GlobalTrackingGeometryRecord>().get(geometryH);
iSetup.get<TrackingComponentsRecord>().get(estimatorName_,estimatorH);
iEvent.getByToken(measurementTrackerTag_, measurementTrackerH);

/// Read L2 track collection
edm::Handle<reco::TrackCollection> l2TrackCol;
iEvent.getByToken(src_, l2TrackCol);

// The product:
std::unique_ptr<std::vector<TrajectorySeed> > result(new std::vector<TrajectorySeed>());

// Get vector of Detector layers once:
std::vector<BarrelDetLayer const*> const& tob = measurementTracker_->geometricSearchTracker()->tobLayers();
std::vector<ForwardDetLayer const*> const& tecPositive = measurementTracker_->geometricSearchTracker()->posTecLayers();
std::vector<ForwardDetLayer const*> const& tecNegative = measurementTracker_->geometricSearchTracker()->negTecLayers();
std::vector<BarrelDetLayer const*> const& tob = measurementTrackerH->geometricSearchTracker()->tobLayers();
std::vector<ForwardDetLayer const*> const& tecPositive = measurementTrackerH->geometricSearchTracker()->posTecLayers();
std::vector<ForwardDetLayer const*> const& tecNegative = measurementTrackerH->geometricSearchTracker()->negTecLayers();
edm::ESHandle<TrackerTopology> tTopo_handle;
iSetup.get<TrackerTopologyRcd>().get(tTopo_handle);
const TrackerTopology* tTopo = tTopo_handle.product();

// Get the suitable propagators:
std::unique_ptr<Propagator> propagatorAlong = SetPropagationDirection(*propagatorAlong_,alongMomentum);
std::unique_ptr<Propagator> propagatorOpposite = SetPropagationDirection(*propagatorOpposite_,oppositeToMomentum);
std::unique_ptr<Propagator> propagatorAlong = SetPropagationDirection(*propagatorAlongH,alongMomentum);
std::unique_ptr<Propagator> propagatorOpposite = SetPropagationDirection(*propagatorOppositeH,oppositeToMomentum);

edm::ESHandle<Propagator> SmartOpposite;
edm::ESHandle<Propagator> SHPOpposite;
Expand All @@ -88,13 +107,13 @@ void TSGForOI::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
std::unique_ptr<std::vector<TrajectorySeed> > out(new std::vector<TrajectorySeed>());
LogTrace("TSGForOI") << "TSGForOI::produce: L2 muon pT, eta, phi --> " << l2->pt() << " , " << l2->eta() << " , " << l2->phi() << endl;

FreeTrajectoryState fts = trajectoryStateTransform::initialFreeState(*l2, magfield_.product());
dummyPlane_->move(fts.position() - dummyPlane_->position());
TrajectoryStateOnSurface tsosAtIP = TrajectoryStateOnSurface(fts, *dummyPlane_);
FreeTrajectoryState fts = trajectoryStateTransform::initialFreeState(*l2, magfieldH.product());
dummyPlane->move(fts.position() - dummyPlane->position());
TrajectoryStateOnSurface tsosAtIP = TrajectoryStateOnSurface(fts, *dummyPlane);
LogTrace("TSGForOI") << "TSGForOI::produce: Created TSOSatIP: " << tsosAtIP << std::endl;

// get the TSOS on the innermost layer of the L2.
TrajectoryStateOnSurface tsosAtMuonSystem = trajectoryStateTransform::innerStateOnSurface(*l2, *geometry_, magfield_.product());
TrajectoryStateOnSurface tsosAtMuonSystem = trajectoryStateTransform::innerStateOnSurface(*l2, *geometryH, magfieldH.product());
LogTrace("TSGForOI") << "TSGForOI::produce: Created TSOSatMuonSystem: " << tsosAtMuonSystem <<endl;

if (useHitLessSeeds_){ //
Expand All @@ -114,43 +133,46 @@ void TSGForOI::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
dist = sqrt(deta*deta+dphi*dphi);
}
if (dist>tsosDiff_){
++numOfMaxSeeds_; // add a hit-based seed
++numOfMaxSeeds; // add a hit-based seed
}
}

numSeedsMade_=0;
analysedL2_ = false;
foundHitlessSeed_ = false;
numSeedsMade=0;
analysedL2 = false;
foundHitlessSeed = false;

// BARREL
if (std::abs(l2->eta()) < maxEtaForTOB_) {
layerCount_ = 0;
layerCount = 0;
for (auto it=tob.rbegin(); it!=tob.rend(); ++it) { //This goes from outermost to innermost layer
LogTrace("TSGForOI") << "TSGForOI::produce: looping in TOB layer " << layerCount_ << endl;
findSeedsOnLayer(tTopo, **it, tsosAtIP, *(propagatorAlong.get()), *(propagatorOpposite.get()), l2, out);
LogTrace("TSGForOI") << "TSGForOI::produce: looping in TOB layer " << layerCount << endl;
findSeedsOnLayer(tTopo, **it, tsosAtIP, *(propagatorAlong.get()), *(propagatorOpposite.get()), l2,
estimatorH, measurementTrackerH, numSeedsMade, numOfMaxSeeds, layerCount, foundHitlessSeed, analysedL2, out);
}
}

// Reset Number of seeds if in overlap region:
if (std::abs(l2->eta())>minEtaForTEC_ && std::abs(l2->eta())<maxEtaForTOB_){
numSeedsMade_=0;
numSeedsMade=0;
}

// ENDCAP+
if (l2->eta() > minEtaForTEC_) {
layerCount_ = 0;
layerCount = 0;
for (auto it=tecPositive.rbegin(); it!=tecPositive.rend(); ++it) {
LogTrace("TSGForOI") << "TSGForOI::produce: looping in TEC+ layer " << layerCount_ << endl;
findSeedsOnLayer(tTopo, **it, tsosAtIP, *(propagatorAlong.get()), *(propagatorOpposite.get()), l2, out);
LogTrace("TSGForOI") << "TSGForOI::produce: looping in TEC+ layer " << layerCount << endl;
findSeedsOnLayer(tTopo, **it, tsosAtIP, *(propagatorAlong.get()), *(propagatorOpposite.get()), l2,
estimatorH, measurementTrackerH, numSeedsMade, numOfMaxSeeds, layerCount, foundHitlessSeed, analysedL2, out);
}
}

// ENDCAP-
if (l2->eta() < -minEtaForTEC_) {
layerCount_ = 0;
layerCount = 0;
for (auto it=tecNegative.rbegin(); it!=tecNegative.rend(); ++it) {
LogTrace("TSGForOI") << "TSGForOI::produce: looping in TEC- layer " << layerCount_ << endl;
findSeedsOnLayer(tTopo, **it, tsosAtIP, *(propagatorAlong.get()), *(propagatorOpposite.get()), l2, out);
LogTrace("TSGForOI") << "TSGForOI::produce: looping in TEC- layer " << layerCount << endl;
findSeedsOnLayer(tTopo, **it, tsosAtIP, *(propagatorAlong.get()), *(propagatorOpposite.get()), l2,
estimatorH, measurementTrackerH, numSeedsMade, numOfMaxSeeds, layerCount, foundHitlessSeed, analysedL2, out);
}
}

Expand All @@ -170,21 +192,28 @@ void TSGForOI::findSeedsOnLayer(
const Propagator& propagatorAlong,
const Propagator& propagatorOpposite,
const reco::TrackRef l2,
std::unique_ptr<std::vector<TrajectorySeed> >& out) {
edm::ESHandle<Chi2MeasurementEstimatorBase>& estimatorH,
edm::Handle<MeasurementTrackerEvent>& measurementTrackerH,
unsigned int& numSeedsMade,
unsigned int& numOfMaxSeeds,
unsigned int& layerCount,
bool& foundHitlessSeed,
bool& analysedL2,
std::unique_ptr<std::vector<TrajectorySeed> >& out) const{

if (numSeedsMade_>numOfMaxSeeds_) return;
LogTrace("TSGForOI") << "TSGForOI::findSeedsOnLayer: numSeedsMade = " << numSeedsMade_ << " , layerCount = " << layerCount_ << endl;
if (numSeedsMade>numOfMaxSeeds) return;
LogTrace("TSGForOI") << "TSGForOI::findSeedsOnLayer: numSeedsMade = " << numSeedsMade << " , layerCount = " << layerCount << endl;

double errorSFHits_=1.0;
double errorSFHitless_=1.0;
if (!adjustErrorsDynamicallyForHits_) errorSFHits_ = fixedErrorRescalingForHits_;
if (!adjustErrorsDynamicallyForHitless_) errorSFHitless_ = fixedErrorRescalingForHitless_;
double errorSFHits=1.0;
double errorSFHitless=1.0;
if (!adjustErrorsDynamicallyForHits_) errorSFHits = fixedErrorRescalingForHits_;
if (!adjustErrorsDynamicallyForHitless_) errorSFHitless = fixedErrorRescalingForHitless_;

// Hitless:
if (useHitLessSeeds_ && !foundHitlessSeed_) {
if (useHitLessSeeds_ && !foundHitlessSeed) {
LogTrace("TSGForOI") << "TSGForOI::findSeedsOnLayer: Start hitless" << endl;
std::vector< GeometricSearchDet::DetWithState > dets;
layer.compatibleDetsV(tsosAtIP, propagatorAlong, *estimator_, dets);
layer.compatibleDetsV(tsosAtIP, propagatorAlong, *estimatorH, dets);
if (dets.size()>0) {
auto const& detOnLayer = dets.front().first;
auto const& tsosOnLayer = dets.front().second;
Expand All @@ -194,31 +223,31 @@ void TSGForOI::findSeedsOnLayer(
}
else{
// calculate SF from L2 (only once -- if needed)
if (!analysedL2_ && adjustErrorsDynamicallyForHitless_) {
errorSFHitless_=calculateSFFromL2(l2);
analysedL2_=true;
if (!analysedL2 && adjustErrorsDynamicallyForHitless_) {
errorSFHitless=calculateSFFromL2(l2);
analysedL2=true;
}

dets.front().second.rescaleError(errorSFHitless_);
dets.front().second.rescaleError(errorSFHitless);
PTrajectoryStateOnDet const& ptsod = trajectoryStateTransform::persistentState(tsosOnLayer,detOnLayer->geographicalId().rawId());
TrajectorySeed::recHitContainer rHC;
out->push_back(TrajectorySeed(ptsod,rHC,oppositeToMomentum));
LogTrace("TSGForOI") << "TSGForOI::findSeedsOnLayer: TSOD (Hitless) done " << endl;
foundHitlessSeed_=true;
foundHitlessSeed=true;
}
numSeedsMade_=out->size();
numSeedsMade=out->size();
}
}
// numSeedsMade_=out->size();
// numSeedsMade=out->size();

// Hits:
if (layerCount_>numOfLayersToTry_) return;
if (layerCount>numOfLayersToTry_) return;
LogTrace("TSGForOI") << "TSGForOI::findSeedsOnLayer: Start Hits" <<endl;
if (makeSeedsFromHits(tTopo, layer, tsosAtIP, *out, propagatorAlong, *measurementTracker_, errorSFHits_)) ++layerCount_;
numSeedsMade_=out->size();
if (makeSeedsFromHits(tTopo, layer, tsosAtIP, *out, propagatorAlong, *measurementTrackerH, estimatorH, errorSFHits)) ++layerCount;
numSeedsMade=out->size();
}

double TSGForOI::calculateSFFromL2(const reco::TrackRef track){
double TSGForOI::calculateSFFromL2(const reco::TrackRef track) const{

double theSF=1.0;
// L2 direction vs pT blowup - as was previously done:
Expand Down Expand Up @@ -250,14 +279,15 @@ int TSGForOI::makeSeedsFromHits(
std::vector<TrajectorySeed> &out,
const Propagator& propagatorAlong,
const MeasurementTrackerEvent &measurementTracker,
const double errorSF) {
edm::ESHandle<Chi2MeasurementEstimatorBase>& estimatorH,
const double errorSF) const{

// Error Rescaling:
TrajectoryStateOnSurface onLayer(tsosAtIP);
onLayer.rescaleError(errorSF);

std::vector< GeometricSearchDet::DetWithState > dets;
layer.compatibleDetsV(onLayer, propagatorAlong, *estimator_, dets);
layer.compatibleDetsV(onLayer, propagatorAlong, *estimatorH, dets);

// Find Measurements on each DetWithState:
LogTrace("TSGForOI") << "TSGForOI::findSeedsOnLayer: find measurements on each detWithState " << dets.size() << endl;
Expand All @@ -269,7 +299,7 @@ int TSGForOI::makeSeedsFromHits(
}
if (!it->second.isValid()) continue; //Skip if TSOS is not valid

std::vector < TrajectoryMeasurement > mymeas = det.fastMeasurements(it->second, onLayer, propagatorAlong, *estimator_); //Second TSOS is not used
std::vector < TrajectoryMeasurement > mymeas = det.fastMeasurements(it->second, onLayer, propagatorAlong, *estimatorH); //Second TSOS is not used
for (std::vector<TrajectoryMeasurement>::const_iterator it2 = mymeas.begin(), ed2 = mymeas.end(); it2 != ed2; ++it2) {
if (it2->recHit()->isValid()) meas.push_back(*it2); //Only save those which are valid
}
Expand Down
64 changes: 29 additions & 35 deletions RecoMuon/TrackerSeedGenerator/plugins/TSGForOI.h
View file
Open in desktop
Expand Up @@ -8,7 +8,7 @@
*/

#include "DataFormats/TrackReco/interface/Track.h"
#include "FWCore/Framework/interface/stream/EDProducer.h"
#include "FWCore/Framework/interface/global/EDProducer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESHandle.h"
Expand All @@ -30,18 +30,18 @@
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "TrackingTools/GeomPropagators/interface/StateOnTrackerBound.h"

class TSGForOI : public edm::stream::EDProducer<> {
class TSGForOI : public edm::global::EDProducer<> {
public:
explicit TSGForOI(const edm::ParameterSet & iConfig);
virtual ~TSGForOI();
static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
virtual void produce(edm::Event& iEvent, const edm::EventSetup& iSetup) override;
virtual void produce(edm::StreamID sid, edm::Event& iEvent, const edm::EventSetup& iSetup) const override;
private:
/// Labels for input collections
const edm::EDGetTokenT<reco::TrackCollection> src_;

/// Maximum number of seeds for each L2
unsigned int numOfMaxSeeds_;
const unsigned int numOfMaxSeedsParam_;

/// How many layers to try
const unsigned int numOfLayersToTry_;
Expand Down Expand Up @@ -73,11 +73,8 @@ class TSGForOI : public edm::stream::EDProducer<> {
/// Switch ON to use Stereo layers instead of using every layer in TEC.
const bool useStereoLayersInTEC_;

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

/// KFUpdator defined in constructor
std::unique_ptr<TrajectoryStateUpdator> updator_;
const std::unique_ptr<TrajectoryStateUpdator> updator_;

const edm::EDGetTokenT<MeasurementTrackerEvent> measurementTrackerTag_;

Expand All @@ -90,41 +87,38 @@ class TSGForOI : public edm::stream::EDProducer<> {
const double tsosDiff_;

/// Counters and flags for the implementation
bool analysedL2_;
bool foundHitlessSeed_;
unsigned int numSeedsMade_;
unsigned int layerCount_;

std::string theCategory;
edm::ESHandle<MagneticField> magfield_;
edm::ESHandle<Propagator> propagatorAlong_;
edm::ESHandle<Propagator> propagatorOpposite_;
edm::ESHandle<GlobalTrackingGeometry> geometry_;
edm::Handle<MeasurementTrackerEvent> measurementTracker_;
edm::ESHandle<Chi2MeasurementEstimatorBase> estimator_;
const std::string theCategory;

/// Function to find seeds on a given layer
void findSeedsOnLayer(
const TrackerTopology* tTopo,
const GeometricSearchDet &layer,
const TrajectoryStateOnSurface &tsosAtIP,
const Propagator& propagatorAlong,
const Propagator& propagatorOpposite,
const reco::TrackRef l2,
std::unique_ptr<std::vector<TrajectorySeed> >& seeds);
const TrackerTopology* tTopo,
const GeometricSearchDet &layer,
const TrajectoryStateOnSurface &tsosAtIP,
const Propagator& propagatorAlong,
const Propagator& propagatorOpposite,
const reco::TrackRef l2,
edm::ESHandle<Chi2MeasurementEstimatorBase>& estimator_,
edm::Handle<MeasurementTrackerEvent>& measurementTrackerH,
unsigned int& numSeedsMade,
unsigned int& numOfMaxSeeds,
unsigned int& layerCount,
bool& foundHitlessSeed,
bool& analysedL2,
std::unique_ptr<std::vector<TrajectorySeed> >& out) const;

/// Function used to calculate the dynamic error SF by analysing the L2
double calculateSFFromL2(const reco::TrackRef track);
double calculateSFFromL2(const reco::TrackRef track) const;

/// Function to find hits on layers and create seeds from updated TSOS
int makeSeedsFromHits(
const TrackerTopology* tTopo,
const GeometricSearchDet &layer,
const TrajectoryStateOnSurface &state,
std::vector<TrajectorySeed> &out,
const Propagator& propagatorAlong,
const MeasurementTrackerEvent &mte,
double errorSF);
const TrackerTopology* tTopo,
const GeometricSearchDet &layer,
const TrajectoryStateOnSurface &tsosAtIP,
std::vector<TrajectorySeed> &out,
const Propagator& propagatorAlong,
const MeasurementTrackerEvent &measurementTracker,
edm::ESHandle<Chi2MeasurementEstimatorBase>& estimator_,
const double errorSF) const;


};
Expand Down