Skip to content

Commit

Permalink
Merge pull request cms-sw#25199 from slava77/CMSSW_10_4_X_2018-11-09-…
Browse files Browse the repository at this point in the history 
…1100/sign1055/eleSeedVector

elide newCombinedSeeds in electron seeding
  • Loading branch information
@cmsbuild
cmsbuild authored Nov 12, 2018
2 parents ef48738 + 2e345f4 commit 524e287
Show file tree
Hide file tree
Showing 7 changed files with 187 additions and 159 deletions.
11 changes: 6 additions & 5 deletions RecoEgamma/EgammaElectronAlgos/interface/ElectronSeedGenerator.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,13 +69,15 @@ class ElectronSeedGenerator
void run(
edm::Event&, const edm::EventSetup& setup,
const reco::SuperClusterRefVector &, const std::vector<float> & hoe1s, const std::vector<float> & hoe2s,
TrajectorySeedCollection *seeds, reco::ElectronSeedCollection&);
const std::vector<const TrajectorySeedCollection *>& seedsV, reco::ElectronSeedCollection&);

private:

void seedsFromThisCluster( edm::Ref<reco::SuperClusterCollection> seedCluster, float hoe1, float hoe2, reco::ElectronSeedCollection & out, const TrackerTopology *tTopo ) ;
void seedsFromRecHits( std::vector<std::pair<RecHitWithDist,ConstRecHitPointer> > & elePixelHits, PropagationDirection & dir, const GlobalPoint & vertexPos, const reco::ElectronSeed::CaloClusterRef & cluster, reco::ElectronSeedCollection & out, bool positron ) ;
void seedsFromTrajectorySeeds( const std::vector<SeedWithInfo> & elePixelSeeds, const reco::ElectronSeed::CaloClusterRef & cluster, float hoe1, float hoe2, reco::ElectronSeedCollection & out, bool positron ) ;
void seedsFromRecHits( std::vector<std::pair<RecHitWithDist,ConstRecHitPointer> > & elePixelHits, PropagationDirection & dir,
const GlobalPoint & vertexPos, const reco::ElectronSeed::CaloClusterRef & cluster, reco::ElectronSeedCollection & out, bool positron ) ;
void seedsFromTrajectorySeeds( const std::vector<SeedWithInfo> & elePixelSeeds, const reco::ElectronSeed::CaloClusterRef & cluster,
float hoe1, float hoe2, reco::ElectronSeedCollection & out, bool positron ) ;
void addSeed( reco::ElectronSeed & seed, const SeedWithInfo * info, bool positron, reco::ElectronSeedCollection & out ) ;
bool prepareElTrackSeed( ConstRecHitPointer outerhit,ConstRecHitPointer innerhit, const GlobalPoint & vertexPos) ;

Expand Down Expand Up @@ -109,8 +111,7 @@ class ElectronSeedGenerator
PixelHitMatcher *myMatchEle;
PixelHitMatcher *myMatchPos;

// edm::Handle<TrajectorySeedCollection> theInitialSeedColl;
TrajectorySeedCollection* theInitialSeedColl;
const std::vector<const TrajectorySeedCollection*>* theInitialSeedCollV = nullptr;

edm::ESHandle<MagneticField> theMagField;
edm::ESHandle<TrackerGeometry> theTrackerGeometry;
Expand Down
2 changes: 1 addition & 1 deletion RecoEgamma/EgammaElectronAlgos/interface/PixelHitMatcher.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -178,7 +178,7 @@ class PixelHitMatcher

std::vector<SeedWithInfo>
compatibleSeeds
( TrajectorySeedCollection * seeds, const GlobalPoint & xmeas,
( const std::vector<const TrajectorySeedCollection *>& seedsV, const GlobalPoint & xmeas,
const GlobalPoint & vprim, float energy, float charge ) ;


Expand Down
15 changes: 6 additions & 9 deletions RecoEgamma/EgammaElectronAlgos/src/ElectronSeedGenerator.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -220,9 +220,9 @@ bool equivalent( const TrajectorySeed & s1, const TrajectorySeed & s2 )
void ElectronSeedGenerator::run
( edm::Event & e, const edm::EventSetup & setup,
const reco::SuperClusterRefVector & sclRefs, const std::vector<float> & hoe1s, const std::vector<float> & hoe2s,
TrajectorySeedCollection * seeds, reco::ElectronSeedCollection & out )
const std::vector<const TrajectorySeedCollection *>& seedsV, reco::ElectronSeedCollection & out )
{
theInitialSeedColl = seeds ;
theInitialSeedCollV = &seedsV;
// bool duplicateTrajectorySeeds =false ;
// unsigned int i,j ;
// for (i=0;i<seeds->size();++i)
Expand Down Expand Up @@ -260,9 +260,6 @@ void ElectronSeedGenerator::run
myMatchEle->setEvent(*data);
myMatchPos->setEvent(*data);

// get initial TrajectorySeeds if necessary
// if (fromTrackerSeeds_) e.getByToken(initialSeeds_, theInitialSeedColl);

// get the beamspot from the Event:
//e.getByType(theBeamSpot);
e.getByToken(beamSpotTag_,theBeamSpot);
Expand Down Expand Up @@ -371,11 +368,11 @@ void ElectronSeedGenerator::seedsFromThisCluster
{
// try electron
std::vector<SeedWithInfo> elePixelSeeds
= myMatchEle->compatibleSeeds(theInitialSeedColl,clusterPos,vertexPos,clusterEnergy,-1.) ;
= myMatchEle->compatibleSeeds(*theInitialSeedCollV,clusterPos,vertexPos,clusterEnergy,-1.) ;
seedsFromTrajectorySeeds(elePixelSeeds,caloCluster,hoe1,hoe2,out,false) ;
// try positron
std::vector<SeedWithInfo> posPixelSeeds
= myMatchPos->compatibleSeeds(theInitialSeedColl,clusterPos,vertexPos,clusterEnergy,1.) ;
= myMatchPos->compatibleSeeds(*theInitialSeedCollV,clusterPos,vertexPos,clusterEnergy,1.) ;
seedsFromTrajectorySeeds(posPixelSeeds,caloCluster,hoe1,hoe2,out,true) ;
}

Expand Down Expand Up @@ -425,11 +422,11 @@ void ElectronSeedGenerator::seedsFromThisCluster
{
// try electron
std::vector<SeedWithInfo> elePixelSeeds
= myMatchEle->compatibleSeeds(theInitialSeedColl,clusterPos,vertexPos,clusterEnergy,-1.) ;
= myMatchEle->compatibleSeeds(*theInitialSeedCollV,clusterPos,vertexPos,clusterEnergy,-1.) ;
seedsFromTrajectorySeeds(elePixelSeeds,caloCluster,hoe1,hoe2,out,false) ;
// try positron
std::vector<SeedWithInfo> posPixelSeeds
= myMatchPos->compatibleSeeds(theInitialSeedColl,clusterPos,vertexPos,clusterEnergy,1.) ;
= myMatchPos->compatibleSeeds(*theInitialSeedCollV,clusterPos,vertexPos,clusterEnergy,1.) ;
seedsFromTrajectorySeeds(posPixelSeeds,caloCluster,hoe1,hoe2,out,true) ;
}
}
Expand Down
215 changes: 108 additions & 107 deletions RecoEgamma/EgammaElectronAlgos/src/PixelHitMatcher.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -97,7 +97,7 @@ float PixelHitMatcher::getVertex()

std::vector<SeedWithInfo>
PixelHitMatcher::compatibleSeeds
( TrajectorySeedCollection * seeds, const GlobalPoint & xmeas,
( const std::vector<const TrajectorySeedCollection *>& seedsV, const GlobalPoint & xmeas,
const GlobalPoint & vprim, float energy, float fcharge )
{
typedef std::unordered_map<const GeomDet *, TrajectoryStateOnSurface> DetTsosAssoc;
Expand All @@ -115,11 +115,11 @@ PixelHitMatcher::compatibleSeeds
TrajectoryStateOnSurface tsos(fts, *bpb(fts.position(), fts.momentum()));

std::vector<SeedWithInfo> result ;

//mapTsos_fast_.clear();
unsigned int allSeedsSize = 0;
for (auto const sc : seedsV) allSeedsSize += sc->size();

mapTsos2_fast_.clear();
// mapTsos_fast_.reserve(seeds->size()) ;
mapTsos2_fast_.reserve(seeds->size()) ;
mapTsos2_fast_.reserve(allSeedsSize) ;

// std::vector<TrajectoryStateOnSurface> vTsos(theTrackerGeometry->dets().size());
// TrajectoryStateOnSurface vTsos[theTrackerGeometry->dets().size()];
Expand All @@ -128,110 +128,111 @@ PixelHitMatcher::compatibleSeeds

int iTsos[ndets];
for ( auto & i : iTsos) i=-1;
std::vector<TrajectoryStateOnSurface> vTsos; vTsos.reserve(seeds->size());

for(const auto& seed : *seeds) {
hit_gp_map_.clear();
if( seed.nHits() > 9 ) {
edm::LogWarning("GsfElectronAlgo|UnexpectedSeed") <<"We cannot deal with seeds having more than 9 hits." ;
continue;
}

const TrajectorySeed::range& hits = seed.recHits();
// cache the global points

for( auto it = hits.first; it != hits.second; ++it ) {
hit_gp_map_.emplace_back(it->globalPosition());
}

//iterate on the hits
auto he = hits.second -1;
for( auto it1 = hits.first; it1 < he; ++it1 ) {
if( !it1->isValid() ) continue;
auto idx1 = std::distance(hits.first,it1);
const DetId id1 = it1->geographicalId();
const GeomDet *geomdet1 = it1->det();

auto ix1 = geomdet1->gdetIndex();

/* VI: this generates regression (other cut is just in phi). in my opinion it is safe and makes sense
auto away = geomdet1->position().basicVector().dot(xmeas.basicVector()) <0;
if (away) continue;
*/

const GlobalPoint& hit1Pos = hit_gp_map_[idx1];
auto dt = hit1Pos.x()*xmeas.x()+hit1Pos.y()*xmeas.y();
if (dt<0) continue;
if (dt<phicut*(xmeas_r*hit1Pos.perp())) continue;

if(iTsos[ix1]<0) {
iTsos[ix1] = vTsos.size();
vTsos.push_back(prop1stLayer->propagate(tsos,geomdet1->surface()));
std::vector<TrajectoryStateOnSurface> vTsos; vTsos.reserve(allSeedsSize);

for (const auto seeds : seedsV) {
for(const auto& seed : *seeds) {
hit_gp_map_.clear();
if( seed.nHits() > 9 ) {
edm::LogWarning("GsfElectronAlgo|UnexpectedSeed") <<"We cannot deal with seeds having more than 9 hits." ;
continue;
}
auto tsos1 = &vTsos[iTsos[ix1]];

if( !tsos1->isValid() ) continue;
std::pair<bool, double> est = ( id1.subdetId() % 2 ?
meas1stBLayer.estimate(vprim, *tsos1, hit1Pos) :
meas1stFLayer.estimate(vprim, *tsos1, hit1Pos) );
if( !est.first ) continue;
EleRelPointPair pp1(hit1Pos,tsos1->globalParameters().position(),vprim);
const math::XYZPoint relHit1Pos(hit1Pos-vprim), relTSOSPos(tsos1->globalParameters().position() - vprim);
const int subDet1 = id1.subdetId();
const float dRz1 = ( id1.subdetId()%2 ? pp1.dZ() : pp1.dPerp() );
const float dPhi1 = pp1.dPhi();
// setup our vertex
double zVertex;
if (!useRecoVertex_) {
// we don't know the z vertex position, get it from linear extrapolation
// compute the z vertex from the cluster point and the found pixel hit
const double pxHit1z = hit1Pos.z();
const double pxHit1x = hit1Pos.x();
const double pxHit1y = hit1Pos.y();
const double r1diff = std::sqrt( (pxHit1x-vprim.x())*(pxHit1x-vprim.x()) +
(pxHit1y-vprim.y())*(pxHit1y-vprim.y()) );
const double r2diff = std::sqrt( (xmeas.x()-pxHit1x)*(xmeas.x()-pxHit1x) +
(xmeas.y()-pxHit1y)*(xmeas.y()-pxHit1y) );
zVertex = pxHit1z - r1diff*(xmeas.z()-pxHit1z)/r2diff;
} else {
// here use rather the reco vertex z position
zVertex = vprim.z();

const TrajectorySeed::range& hits = seed.recHits();
// cache the global points

for( auto it = hits.first; it != hits.second; ++it ) {
hit_gp_map_.emplace_back(it->globalPosition());
}
GlobalPoint vertex(vprim.x(),vprim.y(),zVertex);
FreeTrajectoryState fts2 = FTSFromVertexToPointFactory::get(*theMagField, hit1Pos, vertex, energy, charge) ;
// now find the matching hit
for( auto it2 = it1+1; it2 != hits.second; ++it2 ) {
if( !it2->isValid() ) continue;
auto idx2 = std::distance(hits.first,it2);
const DetId id2 = it2->geographicalId();
const GeomDet *geomdet2 = it2->det();
const std::pair<const GeomDet *,GlobalPoint> det_key(geomdet2,hit1Pos);
const TrajectoryStateOnSurface* tsos2;
auto tsos2_itr = mapTsos2_fast_.find(det_key);
if( tsos2_itr != mapTsos2_fast_.end() ) {
tsos2 = &(tsos2_itr->second);
} else {
auto empl_result =
mapTsos2_fast_.emplace(det_key,prop2ndLayer->propagate(fts2,geomdet2->surface()));
tsos2 = &(empl_result.first->second);
}
if( !tsos2->isValid() ) continue;
const GlobalPoint& hit2Pos = hit_gp_map_[idx2];
std::pair<bool,double> est2 = ( id2.subdetId()%2 ?
meas2ndBLayer.estimate(vertex, *tsos2,hit2Pos) :
meas2ndFLayer.estimate(vertex, *tsos2,hit2Pos) );
if (est2.first) {
EleRelPointPair pp2(hit2Pos,tsos2->globalParameters().position(),vertex) ;
const int subDet2 = id2.subdetId();
const float dRz2 = (subDet2%2==1)?pp2.dZ():pp2.dPerp();
const float dPhi2 = pp2.dPhi();
const unsigned char hitsMask = (1<<idx1)|(1<<idx2);
result.push_back(SeedWithInfo(seed,hitsMask,subDet2,dRz2,dPhi2,subDet1,dRz1,dPhi1)) ;
}
}// inner loop on hits
}// outer loop on hits
}// loop on seeds


//iterate on the hits
auto he = hits.second -1;
for( auto it1 = hits.first; it1 < he; ++it1 ) {
if( !it1->isValid() ) continue;
auto idx1 = std::distance(hits.first,it1);
const DetId id1 = it1->geographicalId();
const GeomDet *geomdet1 = it1->det();

auto ix1 = geomdet1->gdetIndex();

/* VI: this generates regression (other cut is just in phi). in my opinion it is safe and makes sense
auto away = geomdet1->position().basicVector().dot(xmeas.basicVector()) <0;
if (away) continue;
*/

const GlobalPoint& hit1Pos = hit_gp_map_[idx1];
auto dt = hit1Pos.x()*xmeas.x()+hit1Pos.y()*xmeas.y();
if (dt<0) continue;
if (dt<phicut*(xmeas_r*hit1Pos.perp())) continue;

if(iTsos[ix1]<0) {
iTsos[ix1] = vTsos.size();
vTsos.push_back(prop1stLayer->propagate(tsos,geomdet1->surface()));
}
auto tsos1 = &vTsos[iTsos[ix1]];

if( !tsos1->isValid() ) continue;
std::pair<bool, double> est = ( id1.subdetId() % 2 ?
meas1stBLayer.estimate(vprim, *tsos1, hit1Pos) :
meas1stFLayer.estimate(vprim, *tsos1, hit1Pos) );
if( !est.first ) continue;
EleRelPointPair pp1(hit1Pos,tsos1->globalParameters().position(),vprim);
const math::XYZPoint relHit1Pos(hit1Pos-vprim), relTSOSPos(tsos1->globalParameters().position() - vprim);
const int subDet1 = id1.subdetId();
const float dRz1 = ( id1.subdetId()%2 ? pp1.dZ() : pp1.dPerp() );
const float dPhi1 = pp1.dPhi();
// setup our vertex
double zVertex;
if (!useRecoVertex_) {
// we don't know the z vertex position, get it from linear extrapolation
// compute the z vertex from the cluster point and the found pixel hit
const double pxHit1z = hit1Pos.z();
const double pxHit1x = hit1Pos.x();
const double pxHit1y = hit1Pos.y();
const double r1diff = std::sqrt( (pxHit1x-vprim.x())*(pxHit1x-vprim.x()) +
(pxHit1y-vprim.y())*(pxHit1y-vprim.y()) );
const double r2diff = std::sqrt( (xmeas.x()-pxHit1x)*(xmeas.x()-pxHit1x) +
(xmeas.y()-pxHit1y)*(xmeas.y()-pxHit1y) );
zVertex = pxHit1z - r1diff*(xmeas.z()-pxHit1z)/r2diff;
} else {
// here use rather the reco vertex z position
zVertex = vprim.z();
}
GlobalPoint vertex(vprim.x(),vprim.y(),zVertex);
FreeTrajectoryState fts2 = FTSFromVertexToPointFactory::get(*theMagField, hit1Pos, vertex, energy, charge) ;
// now find the matching hit
for( auto it2 = it1+1; it2 != hits.second; ++it2 ) {
if( !it2->isValid() ) continue;
auto idx2 = std::distance(hits.first,it2);
const DetId id2 = it2->geographicalId();
const GeomDet *geomdet2 = it2->det();
const std::pair<const GeomDet *,GlobalPoint> det_key(geomdet2,hit1Pos);
const TrajectoryStateOnSurface* tsos2;
auto tsos2_itr = mapTsos2_fast_.find(det_key);
if( tsos2_itr != mapTsos2_fast_.end() ) {
tsos2 = &(tsos2_itr->second);
} else {
auto empl_result =
mapTsos2_fast_.emplace(det_key,prop2ndLayer->propagate(fts2,geomdet2->surface()));
tsos2 = &(empl_result.first->second);
}
if( !tsos2->isValid() ) continue;
const GlobalPoint& hit2Pos = hit_gp_map_[idx2];
std::pair<bool,double> est2 = ( id2.subdetId()%2 ?
meas2ndBLayer.estimate(vertex, *tsos2,hit2Pos) :
meas2ndFLayer.estimate(vertex, *tsos2,hit2Pos) );
if (est2.first) {
EleRelPointPair pp2(hit2Pos,tsos2->globalParameters().position(),vertex) ;
const int subDet2 = id2.subdetId();
const float dRz2 = (subDet2%2==1)?pp2.dZ():pp2.dPerp();
const float dPhi2 = pp2.dPhi();
const unsigned char hitsMask = (1<<idx1)|(1<<idx2);
result.push_back(SeedWithInfo(seed,hitsMask,subDet2,dRz2,dPhi2,subDet1,dRz1,dPhi1)) ;
}
}// inner loop on hits
}// outer loop on hits
}// loop on seeds
}//loop on vector of seeds
mapTsos2_fast_.clear() ;

return result ;
Expand Down
Loading

0 comments on commit 524e287

Please sign in to comment.