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

avoid redundant computations in StripCPE #14344

Merged
merged 4 commits into from May 7, 2016
Merged

avoid redundant computations in StripCPE #14344

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
3baab0d
compiles, moved loop in cpe
May 2, 2016
98808d8
fix misstype
VinInn May 3, 2016
fefdadc
speed also the other methods
VinInn May 3, 2016
e77be1a
use fast cpe also for 1d hits
VinInn May 3, 2016
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
3 changes: 3 additions & 0 deletions DataFormats/Common/interface/DetSetNew.h
View file
Open in desktop
Expand Up @@ -96,6 +96,9 @@ namespace edmNew {
return edm::Ref<typename HandleT::element_type, typename HandleT::element_type::value_type::value_type>( handle.id(), ci, ci - &(container().front()) );
}

unsigned int makeKeyOf(const_iterator ci) const {
return ci - &(container().front());
Copy link
Contributor

Choose a reason for hiding this comment

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

Why not ci - container().begin();? That would seem to me to be easier for a person to understand.

Copy link
Contributor

Choose a reason for hiding this comment

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

Did you have to use front because const_iterator in this class is a pointer and not a std::vector<...>::const_iterator?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

to make sure I am using pointers and not iterators.
it is consistent to the implementation few lines above

}

private:

Expand Down
9 changes: 8 additions & 1 deletion RecoLocalTracker/ClusterParameterEstimator/interface/StripClusterParameterEstimator.h
View file
Open in desktop
Expand Up @@ -13,6 +13,7 @@
#include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementPoint.h"
#include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementError.h"

#include "CommonTools/Utils/interface/DynArray.h"
#include "FWCore/Utilities/interface/Exception.h"


Expand All @@ -25,9 +26,15 @@
class StripClusterParameterEstimator
{
public:
typedef std::pair<LocalPoint,LocalError> LocalValues;
using LocalValues = std::pair<LocalPoint,LocalError>;
using ALocalValues = DynArray<LocalValues>;
using AClusters = DynArray<SiStripCluster const *>;
typedef std::vector<LocalValues> VLocalValues;

virtual void localParameters(AClusters const & clusters, ALocalValues & retValues, const GeomDetUnit& gd, const LocalTrajectoryParameters & ltp) const {
}


virtual LocalValues localParameters( const SiStripCluster&,const GeomDetUnit&) const {
return std::make_pair(LocalPoint(), LocalError());
}
Expand Down
54 changes: 42 additions & 12 deletions RecoLocalTracker/SiStripRecHitConverter/interface/StripCPE.h
View file
Open in desktop
Expand Up @@ -10,14 +10,15 @@
#include "CondFormats/SiStripObjects/interface/SiStripBackPlaneCorrection.h"
#include "CondFormats/SiStripObjects/interface/SiStripConfObject.h"
#include "CondFormats/SiStripObjects/interface/SiStripLatency.h"
#include <ext/hash_map>

class StripTopology;

class StripCPE : public StripClusterParameterEstimator
{
public:
using StripClusterParameterEstimator::localParameters;

StripClusterParameterEstimator::LocalValues localParameters( const SiStripCluster& cl, const GeomDetUnit&) const;
StripClusterParameterEstimator::LocalValues localParameters( const SiStripCluster& cl, const GeomDetUnit&) const override;

StripCPE( edm::ParameterSet & conf,
const MagneticField&,
Expand All @@ -28,6 +29,45 @@ class StripCPE : public StripClusterParameterEstimator
const SiStripLatency&);
LocalVector driftDirection(const StripGeomDetUnit* det) const;

struct Param {
Param() : topology(nullptr) {}
StripTopology const * topology;
LocalVector drift;
float thickness, invThickness,pitch_rel_err2, maxLength;
int nstrips;
float backplanecorrection;
SiStripDetId::ModuleGeometry moduleGeom;
float coveredStrips(const LocalVector&, const LocalPoint&) const;
};


struct AlgoParam {
Param const & p; const LocalTrajectoryParameters & ltp;
SiStripDetId::SubDetector loc; float afullProjection; float corr;
};


virtual StripClusterParameterEstimator::LocalValues
localParameters( const SiStripCluster& cl, AlgoParam const & ap) const {
return std::make_pair(LocalPoint(), LocalError());
}

AlgoParam getAlgoParam(const GeomDetUnit& det, const LocalTrajectoryParameters & ltp) const {

StripCPE::Param const & p = param(det);
SiStripDetId::SubDetector loc = SiStripDetId( det.geographicalId() ).subDetector();

LocalVector track = ltp.momentum();
track *= -p.thickness/track.z();
Copy link
Contributor

Choose a reason for hiding this comment

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

Is this division safe? Could track.z() be zero? We've just been fixing some bugs caused by an invalid TrajectoryStateOnSurface. There might be more.

Copy link
Contributor Author

Choose a reason for hiding this comment

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

was like this before.
In any case things should be checked in advance, not so deep.
Could please remind me which bug you refers to?

Copy link
Contributor

Choose a reason for hiding this comment

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

@VinInn: I was thinking of #14306, which is not related to the code in this PR, but fixes a bug with a bad TrajectoryStateOnSurface. I have a faint recollection of seeing similar issues previously.

How many paths through the code lead to line 61 above? If you have an idea of how to find them all, we could check that the proper validation is being done in advance. But I would prefer bullet-proof code that never has a possibility to propagate NANs and nonsense values.

Copy link
Contributor Author

@VinInn VinInn May 6, 2016 •
edited

Choose a reason for hiding this comment

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

I think we need to have a more in depth discussion about "FPE" and NaN particularly in the context of vectorization and the future move to vector hardware. We cannot afford to protect each single operation. The current accepted wisdom is to let NaN to propagate and trap it at very high level. Maybe we need to invite an expert and give us a lecture on how one manage this type of issues in HPC.
btw #14306 has nothing to do with NaN or "FPE" is just a trajectory not reaching the target.

Copy link
Contributor

Choose a reason for hiding this comment

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

On 5/6/16 8:44 AM, Vincenzo Innocente wrote:

In RecoLocalTracker/SiStripRecHitConverter/interface/StripCPE.h
#14344 (comment):

  • SiStripDetId::SubDetector loc; float afullProjection; float corr;
  • };
  • virtual StripClusterParameterEstimator::LocalValues
  • localParameters( const SiStripCluster& cl, AlgoParam const & ap) const {
  • return std::make_pair(LocalPoint(), LocalError());
  • }
  • AlgoParam getAlgoParam(const GeomDetUnit& det, const LocalTrajectoryParameters & ltp) const {
  • StripCPE::Param const & p = param(det);
  • SiStripDetId::SubDetector loc = SiStripDetId( det.geographicalId() ).subDetector();
  • LocalVector track = ltp.momentum();
  • track *= -p.thickness/track.z();

I think we need to have a more in depth discussion about "FPE" and NaN
particularly in the context of vectorization and the future move to
vector hardware. We cannot afford to protect each single operation. The
current accepted wisdom is to let NaN to propagate and trap it at very
high level.

"A very high level" better be the output of the module or at worst
output of a sequence of modules
that has consumable output downstream.

The problem is more complex for utility/tools which are used in many places.
A fraction of users of the code may not care about execution speed much.
Maybe for that case it's practical to add the checks to avoid FPE
(template to hand over FPE-safe and unsafe interface?).

Maybe we need to invite an expert and give us a lecture on
how one manage this type of issues in HPC.

Do you have a name in mind?
Maybe send a suggestion by email.

—
You are receiving this because you were mentioned.
Reply to this email directly or view it on GitHub
https://github.com/cms-sw/cmssw/pull/14344/files/e77be1aadab9d37f1baf957e190b179f95a1fc39#r62294072

Copy link
Contributor Author

@VinInn VinInn May 6, 2016 •
edited

Choose a reason for hiding this comment

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

Very High level is before storing in the event.
whatever happen in the meantime is irrelevant (including FPEs)

In any case this code was like this before, so this is not the place to argue about it.

We need to investigate. Somebody from Intel or NVidia or a supercomputing center...
In any case it is something that goes beyond CMS


const float fullProjection = p.coveredStrips( track+p.drift, ltp.position());

auto const corr = - 0.5f*(1.f-p.backplanecorrection) * fullProjection
+ 0.5f*p.coveredStrips(track, ltp.position());

return AlgoParam{p,ltp,loc,std::abs(fullProjection),corr};
}

protected:

const bool peakMode_;
Expand All @@ -38,16 +78,6 @@ class StripCPE : public StripClusterParameterEstimator
std::vector<float> xtalk1;
std::vector<float> xtalk2;

struct Param {
Param() : topology(nullptr) {}
StripTopology const * topology;
LocalVector drift;
float thickness, invThickness,pitch_rel_err2, maxLength;
int nstrips;
float backplanecorrection;
SiStripDetId::ModuleGeometry moduleGeom;
float coveredStrips(const LocalVector&, const LocalPoint&) const;
};
Param const & param(const GeomDetUnit& det) const {
return m_Params[det.index()-m_off];
}
Expand Down
11 changes: 10 additions & 1 deletion RecoLocalTracker/SiStripRecHitConverter/interface/StripCPEfromTrackAngle.h
View file
Open in desktop
Expand Up @@ -29,9 +29,18 @@ class StripCPEfromTrackAngle : public StripCPE
Algo m_algo;

public:
using AlgoParam = StripCPE::AlgoParam;
using AClusters = StripClusterParameterEstimator::AClusters;
using ALocalValues = StripClusterParameterEstimator::ALocalValues;

void localParameters(AClusters const & clusters, ALocalValues & retValues, const GeomDetUnit& gd, const LocalTrajectoryParameters &ltp) const override;

StripClusterParameterEstimator::LocalValues
localParameters( const SiStripCluster& cl, AlgoParam const & ap) const override;


StripClusterParameterEstimator::LocalValues
localParameters( const SiStripCluster&, const GeomDetUnit&, const LocalTrajectoryParameters&) const;
localParameters( const SiStripCluster&, const GeomDetUnit&, const LocalTrajectoryParameters&) const override;

float stripErrorSquared(const unsigned N, const float uProj, const SiStripDetId::SubDetector loc ) const ;
float legacyStripErrorSquared(const unsigned N, const float uProj) const;
Expand Down
106 changes: 80 additions & 26 deletions RecoLocalTracker/SiStripRecHitConverter/src/StripCPEfromTrackAngle.cc
View file
Open in desktop
Expand Up @@ -36,6 +36,7 @@ float StripCPEfromTrackAngle::stripErrorSquared(const unsigned N, const float uP
return uerr*uerr;
}


float StripCPEfromTrackAngle::legacyStripErrorSquared(const unsigned N, const float uProj) const {
if unlikely( (float(N)-uProj) > 3.5f )
return float(N*N)/12.f;
Expand All @@ -48,38 +49,91 @@ float StripCPEfromTrackAngle::legacyStripErrorSquared(const unsigned N, const fl
}
}

StripClusterParameterEstimator::LocalValues
StripCPEfromTrackAngle::localParameters( const SiStripCluster& cluster, const GeomDetUnit& det, const LocalTrajectoryParameters& ltp) const {

StripCPE::Param const & p = param(det);
SiStripDetId ssdid = SiStripDetId( det.geographicalId() );

LocalVector track = ltp.momentum();
track *= -p.thickness/track.z();

const unsigned N = cluster.amplitudes().size();
const float fullProjection = p.coveredStrips( track+p.drift, ltp.position());
float uerr2=0;


void StripCPEfromTrackAngle::localParameters(AClusters const & clusters, ALocalValues & retValues, const GeomDetUnit& det, const LocalTrajectoryParameters & ltp) const {


auto const & par = getAlgoParam(det,ltp);
auto const & p = par.p;
auto loc = par.loc;
auto corr = par.corr;
auto afp = par.afullProjection;

auto fill = [&](unsigned int i, float uerr2) {
const float strip = clusters[i]->barycenter() + corr;
retValues[i].first = p.topology->localPosition(strip, ltp.vector());
retValues[i].second = p.topology->localError(strip, uerr2, ltp.vector());
};


switch (m_algo) {
case Algo::chargeCK :
{
auto dQdx = siStripClusterTools::chargePerCM(cluster, ltp, p.invThickness);
uerr2 = dQdx > maxChgOneMIP ? legacyStripErrorSquared(N,std::abs(fullProjection)) : stripErrorSquared( N, std::abs(fullProjection),ssdid.subDetector() );
}
break;
case Algo::legacy :
uerr2 = legacyStripErrorSquared(N,std::abs(fullProjection));
break;
case Algo::mergeCK :
uerr2 = cluster.isMerged() ? legacyStripErrorSquared(N,std::abs(fullProjection)) : stripErrorSquared( N, std::abs(fullProjection),ssdid.subDetector() );
break;
case Algo::chargeCK :
for (auto i=0U; i< clusters.size(); ++i) {
auto dQdx = siStripClusterTools::chargePerCM(*clusters[i], ltp, p.invThickness);
auto N = clusters[i]->amplitudes().size();
auto uerr2 = dQdx > maxChgOneMIP ? legacyStripErrorSquared(N,afp) : stripErrorSquared( N, afp,loc );
fill(i, uerr2);
}
break;
case Algo::legacy :
for (auto i=0U; i< clusters.size(); ++i) {
auto N = clusters[i]->amplitudes().size();
auto uerr2 = legacyStripErrorSquared(N,afp);
fill(i, uerr2);
}
break;
case Algo::mergeCK :
for (auto i=0U; i< clusters.size(); ++i) {
auto N = clusters[i]->amplitudes().size();
auto uerr2 = clusters[i]->isMerged() ? legacyStripErrorSquared(N,afp) : stripErrorSquared( N,afp,loc );
fill(i, uerr2);
}
break;
}

const float strip = cluster.barycenter() - 0.5f*(1.f-p.backplanecorrection) * fullProjection
+ 0.5f*p.coveredStrips(track, ltp.position());

}

StripClusterParameterEstimator::LocalValues
StripCPEfromTrackAngle::localParameters( const SiStripCluster& cluster, AlgoParam const & par) const {
auto const & p = par.p;
auto const & ltp = par.ltp;
auto loc = par.loc;
auto corr = par.corr;
auto afp = par.afullProjection;

float uerr2=0;

auto N = cluster.amplitudes().size();

switch (m_algo) {
case Algo::chargeCK :
{
auto dQdx = siStripClusterTools::chargePerCM(cluster, ltp, p.invThickness);
uerr2 = dQdx > maxChgOneMIP ? legacyStripErrorSquared(N,afp) : stripErrorSquared( N, afp,loc );
}
break;
case Algo::legacy :
uerr2 = legacyStripErrorSquared(N,afp);
break;
case Algo::mergeCK :
uerr2 = cluster.isMerged() ? legacyStripErrorSquared(N,afp) : stripErrorSquared( N, afp,loc );
break;
}

const float strip = cluster.barycenter() + corr;

return std::make_pair( p.topology->localPosition(strip, ltp.vector()),
p.topology->localError(strip, uerr2, ltp.vector()) );
}



StripClusterParameterEstimator::LocalValues
StripCPEfromTrackAngle::localParameters( const SiStripCluster& cluster, const GeomDetUnit& det, const LocalTrajectoryParameters& ltp) const {

auto const & par = getAlgoParam(det,ltp);
return localParameters(cluster,par);
}