Skip to content

Commit

Permalink
Merge 4e5a45a into b8b3c78
Browse files Browse the repository at this point in the history
  • Loading branch information
felicepantaleo authored May 23, 2019
2 parents b8b3c78 + 4e5a45a commit 3551052
Show file tree
Hide file tree
Showing 4 changed files with 420 additions and 372 deletions.
176 changes: 63 additions & 113 deletions RecoLocalCalo/HGCalRecProducers/interface/HGCalCLUEAlgo.h
Original file line number Diff line number Diff line change
Expand Up @@ -7,19 +7,16 @@

#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/HGCRecHit/interface/HGCRecHitCollections.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
#include "Geometry/CaloTopology/interface/HGCalTopology.h"
#include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"

#include "DataFormats/EgammaReco/interface/BasicCluster.h"
#include "DataFormats/Math/interface/Point3D.h"

#include "RecoLocalCalo/HGCalRecProducers/interface/HGCalLayerTiles.h"

#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
#include "RecoLocalCalo/HGCalRecAlgos/interface/KDTreeLinkerAlgoT.h"

// C/C++ headers
#include <set>
Expand All @@ -36,7 +33,6 @@ class HGCalCLUEAlgo : public HGCalClusteringAlgoBase {
(HGCalClusteringAlgoBase::VerbosityLevel)ps.getUntrackedParameter<unsigned int>("verbosity",3),
reco::CaloCluster::undefined),
thresholdW0_(ps.getParameter<std::vector<double> >("thresholdW0")),
positionDeltaRho_c_(ps.getParameter<std::vector<double> >("positionDeltaRho_c")),
vecDeltas_(ps.getParameter<std::vector<double> >("deltac")),
kappa_(ps.getParameter<double>("kappa")),
ecut_(ps.getParameter<double>("ecut")),
Expand All @@ -48,13 +44,14 @@ class HGCalCLUEAlgo : public HGCalClusteringAlgoBase {
nonAgedNoises_(ps.getParameter<edm::ParameterSet>("noises").getParameter<std::vector<double> >("values")),
noiseMip_(ps.getParameter<edm::ParameterSet>("noiseMip").getParameter<double>("value")),
initialized_(false),
points_(2*(maxlayer+1)),
minpos_(2*(maxlayer+1),{ {0.0f,0.0f} }),
maxpos_(2*(maxlayer+1),{ {0.0f,0.0f} }) {}
cells_(2*(maxlayer+1)),
numberOfClustersPerLayer_(2*(maxlayer+1),0)
{}

~HGCalCLUEAlgo() override {}

void populate(const HGCRecHitCollection &hits) override;

// this is the method that will start the clusterisation (it is possible to invoke this method
// more than once - but make sure it is with different hit collections (or else use reset)

Expand All @@ -63,21 +60,16 @@ class HGCalCLUEAlgo : public HGCalClusteringAlgoBase {
// this is the method to get the cluster collection out
std::vector<reco::BasicCluster> getClusters(bool) override;

// use this if you want to reuse the same cluster object but don't want to accumulate clusters
// (hardly useful?)
void reset() override {
clusters_v_.clear();
layerClustersPerLayer_.clear();
for (auto &it : points_) {
it.clear();
std::vector<KDNode>().swap(it);
}
for (unsigned int i = 0; i < minpos_.size(); i++) {
minpos_[i][0] = 0.;
minpos_[i][1] = 0.;
maxpos_[i][0] = 0.;
maxpos_[i][1] = 0.;
for(auto& cl: numberOfClustersPerLayer_)
{
cl = 0;
}

for(auto& cells : cells_)
cells.clear();

}

Density getDensity() override;
Expand All @@ -90,11 +82,6 @@ class HGCalCLUEAlgo : public HGCalClusteringAlgoBase {
2.9,
2.9
});
iDesc.add<std::vector<double>>("positionDeltaRho_c", {
1.3,
1.3,
1.3
});
iDesc.add<std::vector<double>>("deltac", {
1.3,
1.3,
Expand Down Expand Up @@ -122,7 +109,6 @@ class HGCalCLUEAlgo : public HGCalClusteringAlgoBase {
private:
// To compute the cluster position
std::vector<double> thresholdW0_;
std::vector<double> positionDeltaRho_c_;

// The two parameters used to identify clusters
std::vector<double> vecDeltas_;
Expand All @@ -149,96 +135,60 @@ class HGCalCLUEAlgo : public HGCalClusteringAlgoBase {
// initialization bool
bool initialized_;

struct Hexel {
double x;
double y;
double z;
bool isHalfCell;
double weight;
double fraction;
DetId detid;
double rho;
double delta;
int nearestHigher;
int clusterIndex;
float sigmaNoise;
float thickness;
const hgcal::RecHitTools *tools;

Hexel(const HGCRecHit &hit, DetId id_in, bool isHalf, float sigmaNoise_in, float thickness_in,
const hgcal::RecHitTools *tools_in)
: isHalfCell(isHalf),
weight(0.),
fraction(1.0),
detid(id_in),
rho(0.),
delta(0.),
nearestHigher(-1),
clusterIndex(-1),
sigmaNoise(sigmaNoise_in),
thickness(thickness_in),
tools(tools_in) {
const GlobalPoint position(tools->getPosition(detid));
weight = hit.energy();
x = position.x();
y = position.y();
z = position.z();
float outlierDeltaFactor_ = 2.f;

struct CellsOnLayer {
std::vector<DetId> detid;
std::vector<float> x;
std::vector<float> y;

std::vector<float> weight;
std::vector<float> rho;

std::vector<float> delta;
std::vector<int> nearestHigher;
std::vector<int> clusterIndex;
std::vector<float> sigmaNoise;
std::vector< std::vector <int> > followers;
std::vector<bool> isSeed;

void clear()
{
detid.clear();
x.clear();
y.clear();
weight.clear();
rho.clear();
delta.clear();
nearestHigher.clear();
clusterIndex.clear();
sigmaNoise.clear();
followers.clear();
isSeed.clear();
}
Hexel()
: x(0.),
y(0.),
z(0.),
isHalfCell(false),
weight(0.),
fraction(1.0),
detid(),
rho(0.),
delta(0.),
nearestHigher(-1),
clusterIndex(-1),
sigmaNoise(0.),
thickness(0.),
tools(nullptr) {}
bool operator>(const Hexel &rhs) const { return (rho > rhs.rho); }
};

typedef KDTreeLinkerAlgo<Hexel, 2> KDTree;
typedef KDTreeNodeInfoT<Hexel, 2> KDNode;

std::vector<std::vector<std::vector<KDNode> > > layerClustersPerLayer_;

std::vector<size_t> sort_by_delta(const std::vector<KDNode> &v) const {
std::vector<size_t> idx(v.size());
std::iota(std::begin(idx), std::end(idx), 0);
sort(idx.begin(), idx.end(),
[&v](size_t i1, size_t i2) { return v[i1].data.delta > v[i2].data.delta; });
return idx;
}

std::vector<std::vector<KDNode> > points_; // a vector of vectors of hexels, one for each layer
//@@EM todo: the number of layers should be obtained programmatically - the range is 1-n instead
//of 0-n-1...

std::vector<std::array<float, 2> > minpos_;
std::vector<std::array<float, 2> > maxpos_;

// these functions should be in a helper class.
inline double distance2(const Hexel &pt1, const Hexel &pt2) const { // distance squared
const double dx = pt1.x - pt2.x;
const double dy = pt1.y - pt2.y;

std::vector<CellsOnLayer> cells_;

std::vector<int> numberOfClustersPerLayer_;

inline float distance2(int cell1, int cell2, int layerId) const { // distance squared
const float dx = cells_[layerId].x[cell1] - cells_[layerId].x[cell2];
const float dy = cells_[layerId].y[cell1] - cells_[layerId].y[cell2];
return (dx * dx + dy * dy);
} // distance squaredq
inline double distance(const Hexel &pt1,
const Hexel &pt2) const { // 2-d distance on the layer (x-y)
return std::sqrt(distance2(pt1, pt2));
}

inline float distance(int cell1,
int cell2, int layerId) const { // 2-d distance on the layer (x-y)
return std::sqrt(distance2(cell1, cell2, layerId));
}
double calculateLocalDensity(std::vector<KDNode> &, KDTree &,
const unsigned int) const; // return max density
double calculateDistanceToHigher(std::vector<KDNode> &) const;
int findAndAssignClusters(std::vector<KDNode> &, KDTree &, double, KDTreeBox &,
const unsigned int, std::vector<std::vector<KDNode> > &) const;
math::XYZPoint calculatePosition(const std::vector<KDNode> &) const;
void setDensity(const std::vector<KDNode> &nd);

void prepareDataStructures(const unsigned int layerId);
void calculateLocalDensity(const HGCalLayerTiles& lt, const unsigned int layerId, float delta_c); // return max density
void calculateDistanceToHigher(const HGCalLayerTiles& lt, const unsigned int layerId, float delta_c);
int findAndAssignClusters(const unsigned int layerId, float delta_c);
math::XYZPoint calculatePosition(const std::vector<int> &v, const unsigned int layerId) const;
void setDensity(const unsigned int layerId);
};

#endif
67 changes: 67 additions & 0 deletions RecoLocalCalo/HGCalRecProducers/interface/HGCalLayerTiles.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,67 @@
// Authors: Felice Pantaleo - [email protected]
// Date: 03/2019
// Copyright CERN

#ifndef RecoLocalCalo_HGCalRecAlgos_HGCalLayerTiles
#define RecoLocalCalo_HGCalRecAlgos_HGCalLayerTiles

#include "RecoLocalCalo/HGCalRecProducers/interface/HGCalTilesConstants.h"

#include <vector>
#include <array>
#include <cmath>
#include <algorithm>
#include <cassert>

class HGCalLayerTiles {
public:

void fill(const std::vector<float>& x, const std::vector<float>& y) {
auto cellsSize = x.size();
for (unsigned int i = 0; i < cellsSize; ++i) {
tiles_[getGlobalBin(x[i], y[i])].push_back(i);
}
}

int getXBin(float x) const {
constexpr float xRange = hgcaltilesconstants::maxX - hgcaltilesconstants::minX;
static_assert(xRange >= 0.);
constexpr float r = hgcaltilesconstants::nColumns / xRange;
int xBin = (x - hgcaltilesconstants::minX) * r;
xBin = std::clamp(xBin, 0, hgcaltilesconstants::nColumns);
return xBin;
}

int getYBin(float y) const {
constexpr float yRange = hgcaltilesconstants::maxY - hgcaltilesconstants::minY;
static_assert(yRange >= 0.);
constexpr float r = hgcaltilesconstants::nRows / yRange;
int yBin = (y - hgcaltilesconstants::minY) * r;
yBin = std::clamp(yBin, 0, hgcaltilesconstants::nRows);
return yBin;
}

int getGlobalBin(float x, float y) const { return getXBin(x) + getYBin(y) * hgcaltilesconstants::nColumns; }

int getGlobalBinByBin(int xBin, int yBin) const { return xBin + yBin * hgcaltilesconstants::nColumns; }

std::array<int, 4> searchBox(float xMin, float xMax, float yMin, float yMax) const {
int xBinMin = getXBin(xMin);
int xBinMax = getXBin(xMax);
int yBinMin = getYBin(yMin);
int yBinMax = getYBin(yMax);
return std::array<int, 4>({{xBinMin, xBinMax, yBinMin, yBinMax}});
}

void clear() {
for (auto& t : tiles_)
t.clear();
}

const std::vector<int>& operator[](int globalBinId) const { return tiles_[globalBinId]; }

private:
std::array<std::vector<int>, hgcaltilesconstants::nColumns * hgcaltilesconstants::nRows > tiles_;
};

#endif
28 changes: 28 additions & 0 deletions RecoLocalCalo/HGCalRecProducers/interface/HGCalTilesConstants.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,28 @@
// Authors: Felice Pantaleo - [email protected]
// Date: 03/2019
// Copyright CERN

#ifndef RecoLocalCalo_HGCalRecAlgos_interface_HGCalTilesConstants_h
#define RecoLocalCalo_HGCalRecAlgos_interface_HGCalTilesConstants_h

#include <cstdint>
#include <array>

namespace hgcaltilesconstants {

constexpr int32_t ceil(float num) {
return (static_cast<float>(static_cast<int32_t>(num)) == num) ? static_cast<int32_t>(num)
: static_cast<int32_t>(num) + ((num > 0) ? 1 : 0);
}

constexpr float minX = -285.f;
constexpr float maxX = 285.f;
constexpr float minY = -285.f;
constexpr float maxY = 285.f;
constexpr float tileSize = 5.f;
constexpr int nColumns = hgcaltilesconstants::ceil((maxX - minX) / tileSize);
constexpr int nRows = hgcaltilesconstants::ceil((maxY - minY) / tileSize);

} // namespace hgcaltilesconstants

#endif // RecoLocalCalo_HGCalRecAlgos_interface_HGCalTilesConstants_h
Loading

0 comments on commit 3551052

Please sign in to comment.