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.

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

CLUE3D Bug Fixes #36857

Merged
merged 2 commits into from Feb 8, 2022
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
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
23 changes: 23 additions & 0 deletions RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py
Expand Up @@ -32,3 +32,26 @@
)
)
TICL_FEVT.outputCommands.extend(TICL_RECO.outputCommands)

def customiseHGCalOnlyEventContent(process):
def cleanOutputAndSet(outputModule, ticl_outputCommads):
outputModule.outputCommands = ['drop *_*_*_*']
outputModule.outputCommands.extend(ticl_outputCommads)
outputModule.outputCommands.extend(['keep *_HGCalRecHit_*_*',
'keep *_hgcalLayerClusters_*_*',
'keep CaloParticles_mix_*_*',
'keep SimClusters_mix_*_*',
'keep recoTracks_generalTracks_*_*',
'keep recoTrackExtras_generalTracks_*_*',
'keep SimTracks_g4SimHits_*_*',
'keep SimVertexs_g4SimHits_*_*',
'keep *_layerClusterSimClusterAssociationProducer_*_*',
'keep *_layerClusterCaloParticleAssociationProducer_*_*',
])

if hasattr(process, 'FEVTDEBUGEventContent'):
cleanOutputAndSet(process.FEVTDEBUGEventContent, TICL_FEVT.outputCommands)
if hasattr(process, 'FEVTDEBUGHLToutput'):
cleanOutputAndSet(process.FEVTDEBUGHLToutput, TICL_FEVT.outputCommands)

return process
59 changes: 30 additions & 29 deletions RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc
Expand Up @@ -50,9 +50,9 @@ template <typename TILES>
void PatternRecognitionbyCLUE3D<TILES>::dumpTiles(const TILES &tiles) const {
constexpr int nEtaBin = TILES::constants_type_t::nEtaBins;
constexpr int nPhiBin = TILES::constants_type_t::nPhiBins;
auto lastLayerPerSide = (unsigned int)(rhtools_.lastLayer(false));
unsigned int maxLayer = 2 * lastLayerPerSide - 1;
for (unsigned int layer = 0; layer <= maxLayer; layer++) {
auto lastLayerPerSide = static_cast<int>(rhtools_.lastLayer(false));
int maxLayer = 2 * lastLayerPerSide - 1;
for (int layer = 0; layer <= maxLayer; layer++) {
for (int ieta = 0; ieta < nEtaBin; ieta++) {
auto offset = ieta * nPhiBin;
for (int phi = 0; phi < nPhiBin; phi++) {
Expand Down Expand Up @@ -216,13 +216,13 @@ void PatternRecognitionbyCLUE3D<TILES>::makeTracksters(
clusters_[layer].followers.resize(clusters_[layer].x.size());
}

auto lastLayerPerSide = (unsigned int)(rhtools_.lastLayer(false));
unsigned int maxLayer = 2 * lastLayerPerSide - 1;
auto lastLayerPerSide = static_cast<int>(rhtools_.lastLayer(false));
int maxLayer = 2 * lastLayerPerSide - 1;
std::vector<int> numberOfClustersPerLayer(maxLayer, 0);
for (unsigned int i = 0; i <= maxLayer; i++) {
for (int i = 0; i <= maxLayer; i++) {
calculateLocalDensity(input.tiles, i, layerIdx2layerandSoa);
}
for (unsigned int i = 0; i <= maxLayer; i++) {
for (int i = 0; i <= maxLayer; i++) {
calculateDistanceToHigher(input.tiles, i, layerIdx2layerandSoa);
}

Expand Down Expand Up @@ -437,7 +437,7 @@ void PatternRecognitionbyCLUE3D<TILES>::energyRegressionAndID(const std::vector<

template <typename TILES>
void PatternRecognitionbyCLUE3D<TILES>::calculateLocalDensity(
const TILES &tiles, const unsigned int layerId, const std::vector<std::pair<int, int>> &layerIdx2layerandSoa) {
const TILES &tiles, const int layerId, const std::vector<std::pair<int, int>> &layerIdx2layerandSoa) {
constexpr int nEtaBin = TILES::constants_type_t::nEtaBins;
constexpr int nPhiBin = TILES::constants_type_t::nPhiBins;
auto &clustersOnLayer = clusters_[layerId];
Expand All @@ -454,17 +454,17 @@ void PatternRecognitionbyCLUE3D<TILES>::calculateLocalDensity(

for (unsigned int i = 0; i < numberOfClusters; i++) {
// We need to partition the two sides of the HGCAL detector
auto lastLayerPerSide = static_cast<unsigned int>(rhtools_.lastLayer(false)) - 1;
unsigned int minLayer = 0;
unsigned int maxLayer = 2 * lastLayerPerSide - 1;
auto lastLayerPerSide = static_cast<int>(rhtools_.lastLayer(false));
Copy link
Contributor

Choose a reason for hiding this comment

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

This is not the same as previous lastLayerPerSide, and maxLayer as well as a consequence: please confirm that this is wanted

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Ciao @perrotta that's correct, the -1 was a bug as well.

int minLayer = 0;
int maxLayer = 2 * lastLayerPerSide - 1;
if (layerId < lastLayerPerSide) {
minLayer = std::max(layerId - densitySiblingLayers_, minLayer);
maxLayer = std::min(layerId + densitySiblingLayers_, lastLayerPerSide - 1);
} else {
minLayer = std::max(layerId - densitySiblingLayers_, lastLayerPerSide + 1);
minLayer = std::max(layerId - densitySiblingLayers_, lastLayerPerSide);
maxLayer = std::min(layerId + densitySiblingLayers_, maxLayer);
}
for (unsigned int currentLayer = minLayer; currentLayer <= maxLayer; currentLayer++) {
for (int currentLayer = minLayer; currentLayer <= maxLayer; currentLayer++) {
if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > PatternRecognitionAlgoBaseT<TILES>::Advanced) {
edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "RefLayer: " << layerId << " SoaIDX: " << i;
edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "NextLayer: " << currentLayer;
Expand All @@ -474,8 +474,8 @@ void PatternRecognitionbyCLUE3D<TILES>::calculateLocalDensity(
if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > PatternRecognitionAlgoBaseT<TILES>::Advanced) {
edm::LogVerbatim("PatternRecogntionbyCLUE3D") << "onSameLayer: " << onSameLayer;
}
int etaWindow = onSameLayer ? 2 : 1;
int phiWindow = onSameLayer ? 2 : 1;
const int etaWindow = 2;
const int phiWindow = 2;
Copy link
Contributor

@AdrianoDee AdrianoDee Feb 2, 2022

Choose a reason for hiding this comment

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

Just for my knowledge, I most probably have missed some passages: why this was restricted to 1 for layerClusters on different layers?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Ciao @AdrianoDee, the original motivation for having 2 different cuts, within and outside of the layer, is related to the fact that we are clustering in 3D LayerClusters. Therefore there should be no need to further cluster them within every single layer or, if you want, you would need to search a little farther away to find some other clusters to collect on the same layer. The move to 2 was taken from the CA experience: we notice that, especially at high eta values, a single bin in eta,phi can be quite restrictive. Hence to move to just use one value, 2, to cover them all. When we reach the point of including CLUE3D in the reconstruction, we will need, possibly, some further tuning.

int etaBinMin = std::max(tileOnLayer.etaBin(clustersOnLayer.eta[i]) - etaWindow, 0);
int etaBinMax = std::min(tileOnLayer.etaBin(clustersOnLayer.eta[i]) + etaWindow, nEtaBin);
int phiBinMin = tileOnLayer.phiBin(clustersOnLayer.phi[i]) - phiWindow;
Expand Down Expand Up @@ -556,7 +556,7 @@ void PatternRecognitionbyCLUE3D<TILES>::calculateLocalDensity(

template <typename TILES>
void PatternRecognitionbyCLUE3D<TILES>::calculateDistanceToHigher(
const TILES &tiles, const unsigned int layerId, const std::vector<std::pair<int, int>> &layerIdx2layerandSoa) {
const TILES &tiles, const int layerId, const std::vector<std::pair<int, int>> &layerIdx2layerandSoa) {
constexpr int nEtaBin = TILES::constants_type_t::nEtaBins;
constexpr int nPhiBin = TILES::constants_type_t::nPhiBins;
auto &clustersOnLayer = clusters_[layerId];
Expand All @@ -570,9 +570,9 @@ void PatternRecognitionbyCLUE3D<TILES>::calculateDistanceToHigher(
<< tiles[layerId].etaBin(clustersOnLayer.phi[i]);
}
// We need to partition the two sides of the HGCAL detector
auto lastLayerPerSide = static_cast<unsigned int>(rhtools_.lastLayer(false));
unsigned int minLayer = 0;
unsigned int maxLayer = 2 * lastLayerPerSide - 1;
auto lastLayerPerSide = static_cast<int>(rhtools_.lastLayer(false));
int minLayer = 0;
int maxLayer = 2 * lastLayerPerSide - 1;
if (layerId < lastLayerPerSide) {
minLayer = std::max(layerId - densitySiblingLayers_, minLayer);
maxLayer = std::min(layerId + densitySiblingLayers_, lastLayerPerSide - 1);
Expand All @@ -583,10 +583,10 @@ void PatternRecognitionbyCLUE3D<TILES>::calculateDistanceToHigher(
float maxDelta = std::numeric_limits<float>::max();
float i_delta = maxDelta;
std::pair<int, int> i_nearestHigher(-1, -1);
for (unsigned int currentLayer = minLayer; currentLayer <= maxLayer; currentLayer++) {
for (int currentLayer = minLayer; currentLayer <= maxLayer; currentLayer++) {
const auto &tileOnLayer = tiles[currentLayer];
int etaWindow = 2;
int phiWindow = 2;
int etaWindow = 3;
int phiWindow = 3;
int etaBinMin = std::max(tileOnLayer.etaBin(clustersOnLayer.eta[i]) - etaWindow, 0);
int etaBinMax = std::min(tileOnLayer.etaBin(clustersOnLayer.eta[i]) + etaWindow, nEtaBin);
int phiBinMin = tileOnLayer.phiBin(clustersOnLayer.phi[i]) - phiWindow;
Expand Down Expand Up @@ -626,10 +626,10 @@ void PatternRecognitionbyCLUE3D<TILES>::calculateDistanceToHigher(
// update i_nearestHigher
i_nearestHigher = layerandSoa;
}
}
}
}
}
} // End of loop on clusters
} // End of loop on phi bins
} // End of loop on eta bins
} // End of loop on layers

bool foundNearestHigherInEtaPhiCylinder = (i_delta != maxDelta);
if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > PatternRecognitionAlgoBaseT<TILES>::Advanced) {
Expand All @@ -646,7 +646,7 @@ void PatternRecognitionbyCLUE3D<TILES>::calculateDistanceToHigher(
clustersOnLayer.delta[i] = maxDelta;
clustersOnLayer.nearestHigher[i] = {-1, -1};
}
}
} // End of loop on clusters
}

template <typename TILES>
Expand Down Expand Up @@ -681,7 +681,8 @@ int PatternRecognitionbyCLUE3D<TILES>::findAndAssignTracksters(
<< "Found follower on Layer " << layer << " SOAidx: " << i << " attached to cluster on layer: " << lyrIdx
<< " SOAidx: " << soaIdx;
}
clusters_[lyrIdx].followers[soaIdx].emplace_back(layer, i);
if (lyrIdx >= 0)
clusters_[lyrIdx].followers[soaIdx].emplace_back(layer, i);
} else {
if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > PatternRecognitionAlgoBaseT<TILES>::Advanced) {
edm::LogVerbatim("PatternRecogntionbyCLUE3D")
Expand Down Expand Up @@ -718,7 +719,7 @@ void PatternRecognitionbyCLUE3D<TILES>::fillPSetDescription(edm::ParameterSetDes
iDesc.add<bool>("densityOnSameLayer", false);
iDesc.add<double>("criticalEtaPhiDistance", 0.035);
iDesc.add<double>("outlierMultiplier", 2);
iDesc.add<int>("minNumLayerCluster", 5)->setComment("Not Inclusive");
iDesc.add<int>("minNumLayerCluster", 2)->setComment("Not Inclusive");
Copy link
Contributor

Choose a reason for hiding this comment

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

(As above, I have surely missed something). Why is this going from 5 to 2?

Copy link
Contributor Author

Choose a reason for hiding this comment

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

Ciao @AdrianoDee in this case there's no real reason. The idea, in this case, is that the selection will be done on top of this and will likely use more advanced ML techniques. For this reason, it could be better to allow more input to possibly gather a little more in output. I hope this is clear.

Copy link
Contributor

Choose a reason for hiding this comment

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

I see. Thanks @rovere for all the details (here and above).

iDesc.add<std::string>("eid_input_name", "input");
iDesc.add<std::string>("eid_output_name_energy", "output/regressed_energy");
iDesc.add<std::string>("eid_output_name_id", "output/id_probabilities");
Expand Down
4 changes: 2 additions & 2 deletions RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.h
Expand Up @@ -82,8 +82,8 @@ namespace ticl {
c.shrink_to_fit();
}
}
void calculateLocalDensity(const TILES&, const unsigned int layerId, const std::vector<std::pair<int, int>>&);
void calculateDistanceToHigher(const TILES&, const unsigned int layerId, const std::vector<std::pair<int, int>>&);
void calculateLocalDensity(const TILES&, const int layerId, const std::vector<std::pair<int, int>>&);
void calculateDistanceToHigher(const TILES&, const int layerId, const std::vector<std::pair<int, int>>&);
int findAndAssignTracksters(const TILES&, const std::vector<std::pair<int, int>>&);
void dumpClusters(const std::vector<std::pair<int, int>>& layerIdx2layerandSoa, const int) const;
void dumpTracksters(const std::vector<std::pair<int, int>>& layerIdx2layerandSoa,
Expand Down