Skip to content

Commit

Permalink
Merge pull request #27485 from jingyucms/jingyu-lc-dev-base1100pre3
Browse files Browse the repository at this point in the history 
Layer Clustering for HGCal Scintillator
  • Loading branch information
@cmsbuild
cmsbuild committed Aug 2, 2019
2 parents 740b958 + 223cb8c commit d50e504
Show file tree
Hide file tree
Showing 7 changed files with 341 additions and 77 deletions.
6 changes: 6 additions & 0 deletions RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h
View file
Expand Up @@ -44,6 +44,10 @@ namespace hgcal {

bool isHalfCell(const DetId&) const;

bool isSilicon(const DetId&) const;

bool isOnlySilicon(const unsigned int layer) const;

// 4-vector helper functions using GlobalPoint
float getEta(const GlobalPoint& position, const float& vertex_z = 0.) const;
float getPhi(const GlobalPoint& position) const;
Expand All @@ -62,6 +66,7 @@ namespace hgcal {
unsigned int maxNumberOfWafersPerLayer(bool nose = false) const {
return (nose ? maxNumberOfWafersNose_ : maxNumberOfWafersPerLayer_);
}
inline int getScintMaxIphi() const { return bhMaxIphi_; }
inline int getGeometryType() const { return geometryType_; }
bool maskCell(const DetId& id, int corners = 3) const;

Expand All @@ -70,6 +75,7 @@ namespace hgcal {
unsigned int fhOffset_, bhOffset_, bhLastLayer_, fhLastLayer_;
unsigned int maxNumberOfWafersPerLayer_, maxNumberOfWafersNose_;
int geometryType_;
int bhMaxIphi_;
};
} // namespace hgcal

Expand Down
13 changes: 13 additions & 0 deletions RecoLocalCalo/HGCalRecAlgos/src/RecHitTools.cc
View file
Expand Up @@ -89,6 +89,7 @@ void RecHitTools::getEventSetup(const edm::EventSetup& es) {
bhOffset_ =
fhOffset_ + (geomBH->topology().dddConstants()).firstLayer() - (geomEE->topology().dddConstants()).firstLayer();
bhLastLayer_ = bhOffset_ + (geomBH->topology().dddConstants()).layers(true);
bhMaxIphi_ = geomBH->topology().dddConstants().maxCells(true);
} else {
geometryType_ = 0;
geomEE =
Expand Down Expand Up @@ -413,6 +414,18 @@ bool RecHitTools::isHalfCell(const DetId& id) const {
return ishalf;
}

bool RecHitTools::isSilicon(const DetId& id) const {
bool issilicon = false;
if (id.det() == DetId::HGCalEE || id.det() == DetId::HGCalHSi)
issilicon = true;
return issilicon;
}

bool RecHitTools::isOnlySilicon(const unsigned int layer) const {
bool isonlysilicon = (layer % bhLastLayer_) < bhOffset_;
return isonlysilicon;
}

float RecHitTools::getEta(const GlobalPoint& position, const float& vertex_z) const {
GlobalPoint corrected_position = GlobalPoint(position.x(), position.y(), position.z() - vertex_z);
return corrected_position.eta();
Expand Down
37 changes: 28 additions & 9 deletions RecoLocalCalo/HGCalRecProducers/interface/HGCalCLUEAlgo.h
View file
Expand Up @@ -13,6 +13,7 @@

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

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

Expand Down Expand Up @@ -42,6 +43,7 @@ class HGCalCLUEAlgo : public HGCalClusteringAlgoBase {
fcPerEle_(ps.getParameter<double>("fcPerEle")),
nonAgedNoises_(ps.getParameter<edm::ParameterSet>("noises").getParameter<std::vector<double>>("values")),
noiseMip_(ps.getParameter<edm::ParameterSet>("noiseMip").getParameter<double>("noise_MIP")),
use2x2_(ps.getParameter<bool>("use2x2")),
initialized_(false) {}

~HGCalCLUEAlgo() override {}
Expand Down Expand Up @@ -79,6 +81,7 @@ class HGCalCLUEAlgo : public HGCalClusteringAlgoBase {
1.3,
1.3,
5.0,
0.0315, // for scintillator
});
iDesc.add<bool>("dependSensor", true);
iDesc.add<double>("ecut", 3.0);
Expand All @@ -98,6 +101,7 @@ class HGCalCLUEAlgo : public HGCalClusteringAlgoBase {
iDesc.add<edm::ParameterSetDescription>("doseMap", descNestedNoiseMIP);
descNestedNoiseMIP.add<double>("noise_MIP", 1. / 100.);
iDesc.add<edm::ParameterSetDescription>("noiseMip", descNestedNoiseMIP);
iDesc.add<bool>("use2x2", true); // use 2x2 or 3x3 scenario for scint density calculation
}

/// point in the space
Expand Down Expand Up @@ -129,15 +133,20 @@ class HGCalCLUEAlgo : public HGCalClusteringAlgoBase {
std::vector<std::vector<double>> thresholds_;
std::vector<std::vector<double>> v_sigmaNoise_;

bool use2x2_;

// initialization bool
bool initialized_;

float outlierDeltaFactor_ = 2.f;

struct CellsOnLayer {
std::vector<DetId> detid;
std::vector<bool> isSi;
std::vector<float> x;
std::vector<float> y;
std::vector<float> eta;
std::vector<float> phi;

std::vector<float> weight;
std::vector<float> rho;
Expand All @@ -151,8 +160,11 @@ class HGCalCLUEAlgo : public HGCalClusteringAlgoBase {

void clear() {
detid.clear();
isSi.clear();
x.clear();
y.clear();
eta.clear();
phi.clear();
weight.clear();
rho.clear();
delta.clear();
Expand All @@ -168,22 +180,29 @@ class HGCalCLUEAlgo : public HGCalClusteringAlgoBase {

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);
inline float distance2(int cell1, int cell2, int layerId, bool isEtaPhi) const { // distance squared
if (isEtaPhi) {
const float dphi = reco::deltaPhi(cells_[layerId].phi[cell1], cells_[layerId].phi[cell2]);
const float deta = cells_[layerId].eta[cell1] - cells_[layerId].eta[cell2];
return (deta * deta + dphi * dphi);
} else {
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);
}
}

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));
inline float distance(int cell1, int cell2, int layerId, bool isEtaPhi) const { // 2-d distance on the layer (x-y)
return std::sqrt(distance2(cell1, cell2, layerId, isEtaPhi));
}

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);
float delta_c,
float delta_r); // return max density
void calculateDistanceToHigher(const HGCalLayerTiles& lt, const unsigned int layerId, float delta_c, float delta_r);
int findAndAssignClusters(const unsigned int layerId, float delta_c, float delta_r);
math::XYZPoint calculatePosition(const std::vector<int>& v, const unsigned int layerId) const;
void setDensity(const unsigned int layerId);
};
Expand Down
2 changes: 2 additions & 0 deletions RecoLocalCalo/HGCalRecProducers/interface/HGCalClusteringAlgoBase.h
View file
Expand Up @@ -64,6 +64,7 @@ class HGCalClusteringAlgoBase {
lastLayerEE_ = rhtools_.lastLayerEE();
lastLayerFH_ = rhtools_.lastLayerFH();
firstLayerBH_ = rhtools_.firstLayerBH();
scintMaxIphi_ = rhtools_.getScintMaxIphi();
getEventSetupPerAlgorithm(es);
}
inline void setVerbosity(VerbosityLevel the_verbosity) { verbosity_ = the_verbosity; }
Expand All @@ -75,6 +76,7 @@ class HGCalClusteringAlgoBase {
unsigned int lastLayerEE_;
unsigned int lastLayerFH_;
unsigned int firstLayerBH_;
int scintMaxIphi_;

protected:
// The verbosity level
Expand Down
58 changes: 56 additions & 2 deletions RecoLocalCalo/HGCalRecProducers/interface/HGCalLayerTiles.h
View file
Expand Up @@ -14,10 +14,25 @@

class HGCalLayerTiles {
public:
void fill(const std::vector<float>& x, const std::vector<float>& y) {
void fill(const std::vector<float>& x,
const std::vector<float>& y,
const std::vector<float>& eta,
const std::vector<float>& phi,
const std::vector<bool>& isSi) {
auto cellsSize = x.size();
for (unsigned int i = 0; i < cellsSize; ++i) {
tiles_[getGlobalBin(x[i], y[i])].push_back(i);
if (!isSi[i]) {
tiles_[getGlobalBinEtaPhi(eta[i], phi[i])].push_back(i);
// Copy cells in phi=[-3.15,-3.] to the last bin
if (getPhiBin(phi[i]) == mPiPhiBin) {
tiles_[getGlobalBinEtaPhi(eta[i], phi[i] + 2 * M_PI)].push_back(i);
}
// Copy cells in phi=[3.,3.15] to the first bin
if (getPhiBin(phi[i]) == pPiPhiBin) {
tiles_[getGlobalBinEtaPhi(eta[i], phi[i] - 2 * M_PI)].push_back(i);
}
}
}
}

Expand All @@ -39,10 +54,41 @@ class HGCalLayerTiles {
return yBin;
}

int getEtaBin(float eta) const {
constexpr float etaRange = hgcaltilesconstants::maxEta - hgcaltilesconstants::minEta;
static_assert(etaRange >= 0.);
constexpr float r = hgcaltilesconstants::nColumnsEta / etaRange;
int etaBin = (eta - hgcaltilesconstants::minEta) * r;
etaBin = std::clamp(etaBin, 0, hgcaltilesconstants::nColumnsEta);
return etaBin;
}

int getPhiBin(float phi) const {
constexpr float phiRange = hgcaltilesconstants::maxPhi - hgcaltilesconstants::minPhi;
static_assert(phiRange >= 0.);
constexpr float r = hgcaltilesconstants::nRowsPhi / phiRange;
int phiBin = (phi - hgcaltilesconstants::minPhi) * r;
phiBin = std::clamp(phiBin, 0, hgcaltilesconstants::nRowsPhi);
return phiBin;
}

int mPiPhiBin = getPhiBin(-M_PI);
int pPiPhiBin = getPhiBin(M_PI);

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; }

int getGlobalBinEtaPhi(float eta, float phi) const {
return hgcaltilesconstants::nColumns * hgcaltilesconstants::nRows + getEtaBin(eta) +
getPhiBin(phi) * hgcaltilesconstants::nColumnsEta;
}

int getGlobalBinByBinEtaPhi(int etaBin, int phiBin) const {
return hgcaltilesconstants::nColumns * hgcaltilesconstants::nRows + etaBin +
phiBin * hgcaltilesconstants::nColumnsEta;
}

std::array<int, 4> searchBox(float xMin, float xMax, float yMin, float yMax) const {
int xBinMin = getXBin(xMin);
int xBinMax = getXBin(xMax);
Expand All @@ -51,6 +97,14 @@ class HGCalLayerTiles {
return std::array<int, 4>({{xBinMin, xBinMax, yBinMin, yBinMax}});
}

std::array<int, 4> searchBoxEtaPhi(float etaMin, float etaMax, float phiMin, float phiMax) const {
int etaBinMin = getEtaBin(etaMin);
int etaBinMax = getEtaBin(etaMax);
int phiBinMin = getPhiBin(phiMin);
int phiBinMax = getPhiBin(phiMax);
return std::array<int, 4>({{etaBinMin, etaBinMax, phiBinMin, phiBinMax}});
}

void clear() {
for (auto& t : tiles_)
t.clear();
Expand All @@ -59,7 +113,7 @@ class HGCalLayerTiles {
const std::vector<int>& operator[](int globalBinId) const { return tiles_[globalBinId]; }

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

#endif
13 changes: 11 additions & 2 deletions RecoLocalCalo/HGCalRecProducers/interface/HGCalTilesConstants.h
View file
Expand Up @@ -14,14 +14,23 @@ namespace hgcaltilesconstants {
: static_cast<int32_t>(num) + ((num > 0) ? 1 : 0);
}

constexpr float tileSize = 5.f;
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);
constexpr float tileSizeEtaPhi = 0.15f;
constexpr float minEta = -3.f;
constexpr float maxEta = 3.f;
//To properly construct search box for cells in phi=[-3.15,-3.] and [3.,3.15], cells in phi=[3.,3.15] are copied to the first bin and cells in phi=[-3.15,-3.] are copied to the last bin
constexpr float minPhi = -3.3f;
constexpr float maxPhi = 3.3f;
constexpr int nColumnsEta = hgcaltilesconstants::ceil((maxEta - minEta) / tileSizeEtaPhi);
constexpr int nRowsPhi = hgcaltilesconstants::ceil((maxPhi - minPhi) / tileSizeEtaPhi);
constexpr int nTiles = nColumns * nRows + nColumnsEta * nRowsPhi;

} // namespace hgcaltilesconstants

#endif // RecoLocalCalo_HGCalRecAlgos_interface_HGCalTilesConstants_h
#endif // RecoLocalCalo_HGCalRecAlgos_interface_HGCalTilesConstants_h

0 comments on commit d50e504

Please sign in to comment.