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

Improve performance of Cellular Automaton for regional tracking at HLT #17373

Merged
merged 8 commits into from Feb 14, 2017
Merged

Improve performance of Cellular Automaton for regional tracking at HLT #17373

Show file tree
Hide file tree
Changes from 7 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
e737ca1
code cleanup
JanFSchulte Jan 27, 2017
401fea2
small fixes
JanFSchulte Jan 27, 2017
e049096
a bit more cleanup
JanFSchulte Jan 27, 2017
a5a4cdb
code cleanup
JanFSchulte Jan 31, 2017
f150018
code cleanup II
JanFSchulte Jan 31, 2017
55922e5
last cleanup
JanFSchulte Jan 31, 2017
48d82fe
another round of comments from Matti
JanFSchulte Jan 31, 2017
77f9a41
adressing comments from Slava
JanFSchulte Feb 9, 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
20 changes: 9 additions & 11 deletions RecoPixelVertexing/PixelTriplets/plugins/CAHitNtupletEDProducerT.cc
View file
Open in desktop
Expand Up @@ -86,21 +86,19 @@ void CAHitNtupletEDProducerT<T_Generator>::produce(edm::Event& iEvent, const edm
generator_.initEvent(iEvent, iSetup);

LogDebug("CAHitNtupletEDProducer") << "Creating ntuplets for " << regionDoublets.regionSize() << " regions, and " << regionDoublets.layerPairsSize() << " layer pairs";

typename T_Generator::ResultType ntuplets;
ntuplets.reserve(localRA_.upper());

std::vector<typename T_Generator::ResultType> ntuplets;
ntuplets.resize(regionDoublets.regionSize());
for(auto ntuplet : ntuplets) ntuplet.reserve(localRA_.upper());
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

maybe use auto& or auto&&
to be a bit more obvious that it's not a temp copy that's modified.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I added some printouts and It is clear that this loop is now doing operations on the copies.
Please fix.


generator_.hitNtuplets(regionDoublets, ntuplets, iSetup, seedingLayerHits);
int index = 0;
for(const auto& regionLayerPairs: regionDoublets) {
const TrackingRegion& region = regionLayerPairs.region();
auto seedingHitSetsFiller = seedingHitSets->beginRegion(&region);

LogTrace("CAHitNtupletEDProducer") << " starting region";

generator_.hitNtuplets(regionLayerPairs, ntuplets, iSetup, seedingLayerHits);
LogTrace("CAHitNtupletEDProducer") << " created " << ntuplets.size() << " ntuplets";
fillNtuplets(seedingHitSetsFiller, ntuplets);

ntuplets.clear();
fillNtuplets(seedingHitSetsFiller, ntuplets[index]);
ntuplets[index].clear();
index++;
}
localRA_.update(seedingHitSets->size());

Expand Down
274 changes: 236 additions & 38 deletions RecoPixelVertexing/PixelTriplets/plugins/CAHitQuadrupletGenerator.cc
View file
Open in desktop
Expand Up @@ -99,11 +99,8 @@ void CAHitQuadrupletGenerator::fillDescriptions(edm::ParameterSetDescription& de
void CAHitQuadrupletGenerator::initEvent(const edm::Event& ev, const edm::EventSetup& es) {
if (theComparitor) theComparitor->init(ev, es);
}

namespace {
template <typename T_HitDoublets, typename T_GeneratorOrPairsFunction>
void fillGraph(const SeedingLayerSetsHits& layers, CAGraph& g, T_HitDoublets& hitDoublets,
T_GeneratorOrPairsFunction generatorOrPairsFunction) {
void createGraphStructure(const SeedingLayerSetsHits& layers, CAGraph& g) {
for (unsigned int i = 0; i < layers.size(); i++)
{
for (unsigned int j = 0; j < 4; ++j)
Expand Down Expand Up @@ -131,7 +128,41 @@ namespace {
}

}
}
}
}
void clearGraphStructure(const SeedingLayerSetsHits& layers, CAGraph& g) {
g.theLayerPairs.clear();
for (unsigned int i = 0; i < g.theLayers.size(); i++ ){
g.theLayers[i].theInnerLayers.clear();
g.theLayers[i].theInnerLayerPairs.clear();
g.theLayers[i].theOuterLayers.clear();
g.theLayers[i].theOuterLayerPairs.clear();

}

}
template <typename T_HitDoublets, typename T_GeneratorOrPairsFunction>
void fillGraph(const SeedingLayerSetsHits& layers, CAGraph& g, T_HitDoublets& hitDoublets,
T_GeneratorOrPairsFunction generatorOrPairsFunction) {
for (unsigned int i = 0; i < layers.size(); i++)
{
for (unsigned int j = 0; j < 4; ++j)
{
auto vertexIndex = 0;
auto foundVertex = std::find(g.theLayers.begin(), g.theLayers.end(),
layers[i][j].name());
if (foundVertex == g.theLayers.end())
{
vertexIndex = g.theLayers.size() - 1;
}
else
{
vertexIndex = foundVertex - g.theLayers.begin();
}


if (j > 0)
{

auto innerVertex = std::find(g.theLayers.begin(),
Expand Down Expand Up @@ -178,11 +209,12 @@ void CAHitQuadrupletGenerator::hitQuadruplets(const TrackingRegion& region,

CAGraph g;


std::vector<HitDoublets> hitDoublets;


HitPairGeneratorFromLayerPair thePairGenerator(0, 1, &theLayerCache);

createGraphStructure(layers, g);
fillGraph(layers, g, hitDoublets,
[&](const SeedingLayerSetsHits::SeedingLayer& inner,
const SeedingLayerSetsHits::SeedingLayer& outer,
Expand All @@ -202,9 +234,8 @@ void CAHitQuadrupletGenerator::hitQuadruplets(const TrackingRegion& region,
hitQuadruplets(region, result, hitDoubletsPtr, g, es);
theLayerCache.clear();
}

void CAHitQuadrupletGenerator::hitNtuplets(const IntermediateHitDoublets::RegionLayerSets& regionLayerPairs,
OrderedHitSeeds& result,
void CAHitQuadrupletGenerator::hitNtuplets(const IntermediateHitDoublets& regionDoublets,
std::vector<OrderedHitSeeds>& result,
const edm::EventSetup& es,
const SeedingLayerSetsHits& layers) {
CAGraph g;
Expand All @@ -216,21 +247,194 @@ void CAHitQuadrupletGenerator::hitNtuplets(const IntermediateHitDoublets::Region
SeedingLayerSetsHits::LayerIndex outer) {
return pair.innerLayerIndex() == inner && pair.outerLayerIndex() == outer;
};
fillGraph(layers, g, hitDoublets,
[&](const SeedingLayerSetsHits::SeedingLayer& inner,
const SeedingLayerSetsHits::SeedingLayer& outer,
std::vector<const HitDoublets *>& hitDoublets) {
using namespace std::placeholders;
auto found = std::find_if(regionLayerPairs.begin(), regionLayerPairs.end(),
std::bind(layerPairEqual, _1, inner.index(), outer.index()));
if(found != regionLayerPairs.end()) {
hitDoublets.emplace_back(&(found->doublets()));
return true;
}
return false;
});
edm::ESHandle<TrackerTopology> tTopoHand;
es.get<TrackerTopologyRcd>().get(tTopoHand);
const TrackerTopology *tTopo=tTopoHand.product();

const int numberOfHitsInNtuplet = 4;
std::vector<CACell::CAntuplet> foundQuadruplets;



int index =0;
for(const auto& regionLayerPairs: regionDoublets) {

const TrackingRegion& region = regionLayerPairs.region();
hitDoublets.clear();
foundQuadruplets.clear();
if (index == 0){
createGraphStructure(layers, g);
}
else{
clearGraphStructure(layers, g);
}

fillGraph(layers, g, hitDoublets,
[&](const SeedingLayerSetsHits::SeedingLayer& inner,
const SeedingLayerSetsHits::SeedingLayer& outer,
std::vector<const HitDoublets *>& hitDoublets) {
using namespace std::placeholders;
auto found = std::find_if(regionLayerPairs.begin(), regionLayerPairs.end(),
std::bind(layerPairEqual, _1, inner.index(), outer.index()));
if(found != regionLayerPairs.end()) {
hitDoublets.emplace_back(&(found->doublets()));
return true;
}
return false;
});

CellularAutomaton ca(g);

ca.createAndConnectCells(hitDoublets, region, caThetaCut,
caPhiCut, caHardPtCut);

ca.evolve(numberOfHitsInNtuplet);

ca.findNtuplets(foundQuadruplets, numberOfHitsInNtuplet);


const QuantityDependsPtEval maxChi2Eval = maxChi2.evaluator(es);

// re-used thoughout
std::array<float, 4> bc_r;
std::array<float, 4> bc_z;
std::array<float, 4> bc_errZ2;
std::array<GlobalPoint, 4> gps;
std::array<GlobalError, 4> ges;
std::array<bool, 4> barrels;
bool hasAlreadyPushedACandidate = false;
float selectedChi2 = std::numeric_limits<float>::max();
unsigned int previousfourthLayerId = 0;
std::array<unsigned int, 2> previousCellIds ={{0,0}};
unsigned int previousSideId = 0;
int previousSubDetId = 0;

unsigned int numberOfFoundQuadruplets = foundQuadruplets.size();

// Loop over quadruplets
for (unsigned int quadId = 0; quadId < numberOfFoundQuadruplets; ++quadId)
{

auto isBarrel = [](const unsigned id) -> bool
{
return id == PixelSubdetector::PixelBarrel;
};
for(unsigned int i = 0; i< 3; ++i)
{
auto const& ahit = foundQuadruplets[quadId][i]->getInnerHit();
gps[i] = ahit->globalPosition();
ges[i] = ahit->globalPositionError();
barrels[i] = isBarrel(ahit->geographicalId().subdetId());
}

auto const& ahit = foundQuadruplets[quadId][2]->getOuterHit();
gps[3] = ahit->globalPosition();
ges[3] = ahit->globalPositionError();
barrels[3] = isBarrel(ahit->geographicalId().subdetId());
if(caOnlyOneLastHitPerLayerFilter)
{
const auto fourthLayerId = tTopo->layer(ahit->geographicalId());
const auto sideId = tTopo->side(ahit->geographicalId());
const auto subDetId = ahit->geographicalId().subdetId();
const auto isTheSameTriplet = (quadId != 0) && (foundQuadruplets[quadId][0]->getCellId() == previousCellIds[0]) && (foundQuadruplets[quadId][1]->getCellId() == previousCellIds[1]);
const auto isTheSameFourthLayer = (quadId != 0) && (fourthLayerId == previousfourthLayerId) && (subDetId == previousSubDetId) && (sideId == previousSideId);

previousCellIds = {{foundQuadruplets[quadId][0]->getCellId(), foundQuadruplets[quadId][1]->getCellId()}};
previousfourthLayerId = fourthLayerId;


if(!(isTheSameTriplet && isTheSameFourthLayer ))
{
selectedChi2 = std::numeric_limits<float>::max();
hasAlreadyPushedACandidate = false;
}

}
// TODO:
// - 'line' is not used for anything
// - if we decide to always do the circle fit for 4 hits, we don't
// need ThirdHitPredictionFromCircle for the curvature; then we
// could remove extraHitRPhitolerance configuration parameter
PixelRecoLineRZ line(gps[0], gps[2]);
ThirdHitPredictionFromCircle predictionRPhi(gps[0], gps[2], extraHitRPhitolerance);
const float curvature = predictionRPhi.curvature(ThirdHitPredictionFromCircle::Vector2D(gps[1].x(), gps[1].y()));
const float abscurv = std::abs(curvature);
const float thisMaxChi2 = maxChi2Eval.value(abscurv);
if (theComparitor)
{
SeedingHitSet tmpTriplet(foundQuadruplets[quadId][0]->getInnerHit(), foundQuadruplets[quadId][2]->getInnerHit(), foundQuadruplets[quadId][2]->getOuterHit());


if (!theComparitor->compatible(tmpTriplet) )
{
continue;
}
}

float chi2 = std::numeric_limits<float>::quiet_NaN();
// TODO: Do we have any use case to not use bending correction?
if (useBendingCorrection)
{
// Following PixelFitterByConformalMappingAndLine
const float simpleCot = ( gps.back().z() - gps.front().z() ) / (gps.back().perp() - gps.front().perp() );
const float pt = 1.f / PixelRecoUtilities::inversePt(abscurv, es);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

some indentation to offset the code in the if (useBendingCorrection) scope would be useful for readability

for (int i=0; i < 4; ++i)
{
const GlobalPoint & point = gps[i];
const GlobalError & error = ges[i];
bc_r[i] = sqrt( sqr(point.x() - region.origin().x()) + sqr(point.y() - region.origin().y()) );
bc_r[i] += pixelrecoutilities::LongitudinalBendingCorrection(pt, es)(bc_r[i]);
bc_z[i] = point.z() - region.origin().z();
bc_errZ2[i] = (barrels[i]) ? error.czz() : error.rerr(point)*sqr(simpleCot);
}
RZLine rzLine(bc_r, bc_z, bc_errZ2, RZLine::ErrZ2_tag());
chi2 = rzLine.chi2();
}
else
{
RZLine rzLine(gps, ges, barrels);
chi2 = rzLine.chi2();
}
if (edm::isNotFinite(chi2) || chi2 > thisMaxChi2)
{
continue;

}
// TODO: Do we have any use case to not use circle fit? Maybe
// HLT where low-pT inefficiency is not a problem?
if (fitFastCircle)
{
FastCircleFit c(gps, ges);
chi2 += c.chi2();
if (edm::isNotFinite(chi2))
continue;
if (fitFastCircleChi2Cut && chi2 > thisMaxChi2)
continue;
}
if(caOnlyOneLastHitPerLayerFilter)
{
if (chi2 < selectedChi2)
{
selectedChi2 = chi2;

if(hasAlreadyPushedACandidate)
{
result[index].pop_back();

}
result[index].emplace_back(foundQuadruplets[quadId][0]->getInnerHit(), foundQuadruplets[quadId][1]->getInnerHit(),foundQuadruplets[quadId][2]->getInnerHit(), foundQuadruplets[quadId][2]->getOuterHit());
hasAlreadyPushedACandidate = true;

}
}
else
{
result[index].emplace_back(foundQuadruplets[quadId][0]->getInnerHit(), foundQuadruplets[quadId][1]->getInnerHit(), foundQuadruplets[quadId][2]->getInnerHit(), foundQuadruplets[quadId][2]->getOuterHit());
}
}
index++;
}

hitQuadruplets(regionLayerPairs.region(), result, hitDoublets, g, es);
}

void CAHitQuadrupletGenerator::hitQuadruplets(const TrackingRegion& region,
Expand All @@ -257,25 +461,19 @@ void CAHitQuadrupletGenerator::hitQuadruplets(const TrackingRegion& region,

const QuantityDependsPtEval maxChi2Eval = maxChi2.evaluator(es);

// re-used thoughout, need to be vectors because of RZLine interface
// re-used thoughout
std::array<float, 4> bc_r;
std::array<float, 4> bc_z;
std::array<float, 4> bc_errZ2;
std::array<GlobalPoint, 4> gps;
std::array<GlobalError, 4> ges;
std::array<bool, 4> barrels;
unsigned int fourthLayerId = 0;
unsigned int previousfourthLayerId = 0;
int subDetId = 0;
int previousSubDetId = 0;
unsigned int sideId = 0;
unsigned int previousSideId = 0;
std::array<unsigned int, 2> previousCellIds ={{0,0}};
bool isTheSameTriplet = false;
bool isTheSameFourthLayer = false;
bool hasAlreadyPushedACandidate = false;
float selectedChi2 = std::numeric_limits<float>::max();

unsigned int previousfourthLayerId = 0;
std::array<unsigned int, 2> previousCellIds ={{0,0}};
unsigned int previousSideId = 0;
int previousSubDetId = 0;

unsigned int numberOfFoundQuadruplets = foundQuadruplets.size();

Expand All @@ -302,11 +500,11 @@ void CAHitQuadrupletGenerator::hitQuadruplets(const TrackingRegion& region,

if(caOnlyOneLastHitPerLayerFilter)
{
fourthLayerId = tTopo->layer(ahit->geographicalId());
sideId = tTopo->side(ahit->geographicalId());
subDetId = ahit->geographicalId().subdetId();
isTheSameTriplet = (quadId != 0) && (foundQuadruplets[quadId][0]->getCellId() == previousCellIds[0]) && (foundQuadruplets[quadId][1]->getCellId() == previousCellIds[1]);
isTheSameFourthLayer = (quadId != 0) && (fourthLayerId == previousfourthLayerId) && (subDetId == previousSubDetId) && (sideId == previousSideId);
const auto fourthLayerId = tTopo->layer(ahit->geographicalId());
const auto sideId = tTopo->side(ahit->geographicalId());
const auto subDetId = ahit->geographicalId().subdetId();
const auto isTheSameTriplet = (quadId != 0) && (foundQuadruplets[quadId][0]->getCellId() == previousCellIds[0]) && (foundQuadruplets[quadId][1]->getCellId() == previousCellIds[1]);
const auto isTheSameFourthLayer = (quadId != 0) && (fourthLayerId == previousfourthLayerId) && (subDetId == previousSubDetId) && (sideId == previousSideId);

previousCellIds = {{foundQuadruplets[quadId][0]->getCellId(), foundQuadruplets[quadId][1]->getCellId()}};
previousfourthLayerId = fourthLayerId;
Expand Down
4 changes: 2 additions & 2 deletions RecoPixelVertexing/PixelTriplets/plugins/CAHitQuadrupletGenerator.h
View file
Open in desktop
Expand Up @@ -54,8 +54,8 @@ class CAHitQuadrupletGenerator : public HitQuadrupletGenerator {
const edm::Event & ev, const edm::EventSetup& es);

// new-style
void hitNtuplets(const IntermediateHitDoublets::RegionLayerSets& regionLayerPairs,
OrderedHitSeeds& result,
void hitNtuplets(const IntermediateHitDoublets& regionDoublets,
std::vector<OrderedHitSeeds>& result,
const edm::EventSetup& es,
const SeedingLayerSetsHits& layers);

Expand Down