Skip to content

Commit

Permalink
Electron superclustering in TICL using a DNN
Browse files Browse the repository at this point in the history
  • Loading branch information
tcuisset committed Apr 23, 2024
1 parent 50186b3 commit ac56f5d
Show file tree
Hide file tree
Showing 28 changed files with 1,998 additions and 1,547 deletions.
6 changes: 6 additions & 0 deletions CommonTools/RecoAlgos/interface/MultiVectorManager.h
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,12 @@ class MultiVectorManager {

class Iterator {
public:
using iterator_category = std::forward_iterator_tag;
using difference_type = std::ptrdiff_t;
using value_type = T;
using pointer = T*;
using reference = T&;

Iterator(const MultiVectorManager& manager, size_t index) : manager(manager), currentIndex(index) {}

bool operator!=(const Iterator& other) const { return currentIndex != other.currentIndex; }
Expand Down
6 changes: 6 additions & 0 deletions DQMOffline/EGamma/python/electronAnalyzerSequence_cff.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@

import FWCore.ParameterSet.Config as cms

from Configuration.ProcessModifiers.ticl_v5_cff import ticl_v5

mergedSuperClusters = cms.EDProducer("SuperClusterMerger",
src = cms.VInputTag(
# cms.InputTag("correctedHybridSuperClusters"),
Expand All @@ -9,6 +11,10 @@
cms.InputTag("particleFlowSuperClusterECAL","particleFlowSuperClusterECALEndcapWithPreshower")
)
)
ticl_v5.toModify(mergedSuperClusters, src=cms.VInputTag(
cms.InputTag("particleFlowSuperClusterECAL","particleFlowSuperClusterECALBarrel"),
cms.InputTag("ticlEGammaSuperClusterProducer")
))

from DQMOffline.EGamma.electronGeneralAnalyzer_cfi import *
dqmElectronGeneralAnalysis.OutputFolderName = cms.string("Egamma/Electrons/Ele1_General") ;
Expand Down
8 changes: 4 additions & 4 deletions DataFormats/HGCalReco/interface/Trackster.h
Original file line number Diff line number Diff line change
Expand Up @@ -95,10 +95,10 @@ namespace ticl {
removeDuplicates();
zeroProbabilities();
}
inline void fillPCAVariables(Eigen::Vector3d &eigenvalues,
Eigen::Matrix3d &eigenvectors,
Eigen::Vector3d &sigmas,
Eigen::Vector3d &sigmasEigen,
inline void fillPCAVariables(Eigen::Vector3f const &eigenvalues,
Eigen::Matrix3f const &eigenvectors,
Eigen::Vector3f const &sigmas,
Eigen::Vector3f const &sigmasEigen,
size_t pcadimension,
PCAOrdering order) {
int original_index = 0;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,12 @@
allowHGCal = True,
)

from Configuration.ProcessModifiers.ticl_v5_cff import ticl_v5
ticl_v5.toModify(
ecalDrivenElectronSeeds,
endcapSuperClusters = cms.InputTag("ticlEGammaSuperClusterProducer")
)

from Configuration.ProcessModifiers.egamma_lowPt_exclusive_cff import egamma_lowPt_exclusive
egamma_lowPt_exclusive.toModify(ecalDrivenElectronSeeds,
LowPtThreshold =1.0,
Expand Down
5 changes: 5 additions & 0 deletions RecoEgamma/EgammaPhotonProducers/python/photonCore_cfi.py
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,11 @@
pixelSeedProducer = 'electronMergedSeeds',
endcapOnly = True,
)
from Configuration.ProcessModifiers.ticl_v5_cff import ticl_v5
ticl_v5.toModify(
photonCoreHGC,
scIslandEndcapProducer = cms.InputTag("ticlEGammaSuperClusterProducer")
)

islandPhotonCore = photonCore.clone(
scHybridBarrelProducer = "correctedIslandBarrelSuperClusters",
Expand Down
5 changes: 5 additions & 0 deletions RecoHGCal/Configuration/python/RecoHGCal_EventContent_cff.py
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@
'keep *_ticlTrackstersCLUE3DEM_*_*',
'keep *_ticlTrackstersCLUE3DHAD_*_*',
'keep *_ticlTracksterLinks_*_*',
'keep *_ticlTracksterLinksSuperclustering_*_*',
'keep *_ticlCandidate_*_*',
]
)
Expand All @@ -60,6 +61,10 @@
'keep *_ticlSimTracksters_*_*',
'keep *_ticlSimTICLCandidates_*_*',
'keep *_ticlSimTrackstersFromCP_*_*',
'keep *_tracksterSimTracksterAssociationLinkingbyCLUE3DEM_*_*',
'keep *_tracksterSimTracksterAssociationPRbyCLUE3DEM_*_*',
'keep *_tracksterSimTracksterAssociationLinkingSuperclustering_*_*',
'keep *_tracksterSimTracksterAssociationPRSuperclustering_*_*',
)
)

Expand Down
93 changes: 93 additions & 0 deletions RecoHGCal/TICL/interface/SuperclusteringDNNInputs.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,93 @@
/** Produce inputs for the superclustering DNN */
// Author: Theo Cuisset - [email protected]
// Date: 11/2023

#ifndef __RecoHGCal_TICL_SuperclusteringDNNInputs_H__
#define __RecoHGCal_TICL_SuperclusteringDNNInputs_H__

#include <vector>
#include <string>
#include <memory>

namespace ticl {
class Trackster;

// Abstract base class for DNN input preparation.
class AbstractDNNInput {
public:
virtual ~AbstractDNNInput() = default;

virtual unsigned int featureCount() const { return featureNames().size(); };

/** Get name of features. Used for SuperclusteringSampleDumper branch names (inference does not use the names, only the indices)
* The default implementation is meant to be overriden by inheriting classes
*/
virtual std::vector<std::string> featureNames() const {
std::vector<std::string> defaultNames;
defaultNames.reserve(featureCount());
for (unsigned int i = 1; i <= featureCount(); i++) {
defaultNames.push_back(std::string("nb_") + std::to_string(i));
}
return defaultNames;
}

/** Compute feature for seed and candidate pair */
virtual std::vector<float> computeVector(ticl::Trackster const& ts_base, ticl::Trackster const& ts_toCluster) = 0;
};

/* First version of DNN by Alessandro Tarabini. Meant as a DNN equivalent of Moustache algorithm (superclustering algo in ECAL)
Uses features : ['DeltaEta', 'DeltaPhi', 'multi_en', 'multi_eta', 'multi_pt', 'seedEta','seedPhi','seedEn', 'seedPt']
*/
class DNNInputAlessandroV1 : public AbstractDNNInput {
public:
unsigned int featureCount() const override { return 9; }

std::vector<float> computeVector(ticl::Trackster const& ts_base, ticl::Trackster const& ts_toCluster) override;

std::vector<std::string> featureNames() const override {
return {"DeltaEtaBaryc",
"DeltaPhiBaryc",
"multi_en",
"multi_eta",
"multi_pt",
"seedEta",
"seedPhi",
"seedEn",
"seedPt"};
}
};

/* Second version of DNN by Alessandro Tarabini.
Uses features : ['DeltaEta', 'DeltaPhi', 'multi_en', 'multi_eta', 'multi_pt', 'seedEta','seedPhi','seedEn', 'seedPt', 'theta', 'theta_xz_seedFrame', 'theta_yz_seedFrame', 'theta_xy_cmsFrame', 'theta_yz_cmsFrame', 'theta_xz_cmsFrame', 'explVar', 'explVarRatio']
*/
class DNNInputAlessandroV2 : public AbstractDNNInput {
public:
unsigned int featureCount() const override { return 17; }

std::vector<float> computeVector(ticl::Trackster const& ts_base, ticl::Trackster const& ts_toCluster) override;

std::vector<std::string> featureNames() const override {
return {"DeltaEtaBaryc",
"DeltaPhiBaryc",
"multi_en",
"multi_eta",
"multi_pt",
"seedEta",
"seedPhi",
"seedEn",
"seedPt",
"theta",
"theta_xz_seedFrame",
"theta_yz_seedFrame",
"theta_xy_cmsFrame",
"theta_yz_cmsFrame",
"theta_xz_cmsFrame",
"explVar",
"explVarRatio"};
}
};

std::unique_ptr<AbstractDNNInput> makeDNNInputFromString(std::string dnnVersion);
} // namespace ticl

#endif
16 changes: 14 additions & 2 deletions RecoHGCal/TICL/interface/TracksterLinkingAlgoBase.h
Original file line number Diff line number Diff line change
Expand Up @@ -30,11 +30,22 @@ namespace edm {
class EventSetup;
} // namespace edm

namespace cms {
namespace Ort {
class ONNXRuntime;
}
} // namespace cms

namespace ticl {
class TracksterLinkingAlgoBase {
public:
TracksterLinkingAlgoBase(const edm::ParameterSet& conf, edm::ConsumesCollector)
: algo_verbosity_(conf.getParameter<int>("algo_verbosity")) {}
/** \param conf the configuration of the plugin
* \param onnxRuntime the ONNXRuntime, if onnxModelPath was provided in plugin configuration (nullptr otherwise)
*/
TracksterLinkingAlgoBase(const edm::ParameterSet& conf,
edm::ConsumesCollector,
cms::Ort::ONNXRuntime const* onnxRuntime = nullptr)
: algo_verbosity_(conf.getParameter<int>("algo_verbosity")), onnxRuntime_(onnxRuntime) {}
virtual ~TracksterLinkingAlgoBase(){};

struct Inputs {
Expand Down Expand Up @@ -66,6 +77,7 @@ namespace ticl {

protected:
int algo_verbosity_;
cms::Ort::ONNXRuntime const* onnxRuntime_;
};
} // namespace ticl

Expand Down
1 change: 1 addition & 0 deletions RecoHGCal/TICL/plugins/BuildFile.xml
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,7 @@
<use name="MagneticField/Records"/>
<use name="PhysicsTools/TensorFlow"/>
<use name="PhysicsTools/UtilAlgos"/>
<use name="PhysicsTools/ONNXRuntime" />
<use name="RecoLocalCalo/HGCalRecAlgos"/>
<use name="RecoLocalCalo/HGCalRecProducers"/>
<use name="RecoHGCal/TICL"/>
Expand Down
124 changes: 124 additions & 0 deletions RecoHGCal/TICL/plugins/EGammaSuperclusterProducer.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,124 @@
// Authors : Theo Cuisset <[email protected]>, Shamik Ghosh <[email protected]>
// Date : 01/2024
/* Translates TICL superclusters to ECAL supercluster dataformats (reco::SuperCluster and reco::CaloCluster). Performs similar task as RecoEcal/EgammaCLusterProducers/PFECALSuperClusterProducer */

#include "FWCore/Framework/interface/stream/EDProducer.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ConsumesCollector.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/Utilities/interface/EDGetToken.h"

#include "DataFormats/EgammaReco/interface/SuperClusterFwd.h"
#include "DataFormats/EgammaReco/interface/SuperCluster.h"
#include "DataFormats/EgammaReco/interface/ElectronSeed.h"
#include "DataFormats/EgammaReco/interface/ElectronSeedFwd.h"
#include "DataFormats/HGCalReco/interface/Trackster.h"

class EGammaSuperclusterProducer : public edm::stream::EDProducer<> {
public:
EGammaSuperclusterProducer(const edm::ParameterSet&);

void produce(edm::Event&, const edm::EventSetup&) override;

static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);

private:
edm::EDGetTokenT<ticl::TracksterCollection> ticlSuperClustersToken_;
edm::EDGetTokenT<std::vector<std::vector<unsigned int>>> superClusterLinksToken_;
edm::EDGetTokenT<ticl::TracksterCollection> ticlTrackstersEMToken_;
edm::EDGetTokenT<reco::CaloClusterCollection> layerClustersToken_;
};

EGammaSuperclusterProducer::EGammaSuperclusterProducer(const edm::ParameterSet& ps)
: ticlSuperClustersToken_(consumes<ticl::TracksterCollection>(ps.getParameter<edm::InputTag>("ticlSuperClusters"))),
superClusterLinksToken_(consumes<std::vector<std::vector<unsigned int>>>(
edm::InputTag(ps.getParameter<edm::InputTag>("ticlSuperClusters").label(),
"linkedTracksterIdToInputTracksterId",
ps.getParameter<edm::InputTag>("ticlSuperClusters").process()))),
ticlTrackstersEMToken_(consumes<ticl::TracksterCollection>(ps.getParameter<edm::InputTag>("ticlTrackstersEM"))),
layerClustersToken_(consumes<reco::CaloClusterCollection>(ps.getParameter<edm::InputTag>("layerClusters"))) {
produces<reco::SuperClusterCollection>();
produces<reco::CaloClusterCollection>(); // The CaloCluster corresponding to each EM trackster
}

void EGammaSuperclusterProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
auto const& ticlSuperclusters = iEvent.get(ticlSuperClustersToken_);
auto const& ticlSuperclusterLinks = iEvent.get(superClusterLinksToken_);
auto emTracksters_h = iEvent.getHandle(ticlTrackstersEMToken_);
auto const& emTracksters = *emTracksters_h;
auto const& layerClusters = iEvent.get(layerClustersToken_);
// Output collections :
auto egammaSuperclusters = std::make_unique<reco::SuperClusterCollection>();
auto caloClustersEM = std::make_unique<reco::CaloClusterCollection>();

// Fill reco::CaloCluster collection (1-1 mapping to TICL EM tracksters)
for (ticl::Trackster const& emTrackster : emTracksters) {
std::vector<std::pair<DetId, float>> hitsAndFractions;
int iLC = 0;
std::for_each(std::begin(emTrackster.vertices()), std::end(emTrackster.vertices()), [&](unsigned int lcId) {
const auto fraction = 1.f / emTrackster.vertex_multiplicity(iLC++);
for (const auto& cell : layerClusters[lcId].hitsAndFractions()) {
hitsAndFractions.emplace_back(cell.first, cell.second * fraction);
}
});

reco::CaloCluster& caloCluster = caloClustersEM->emplace_back(
emTrackster.regressed_energy(), // energy
math::XYZPoint(emTrackster.barycenter()), // position
reco::CaloID(reco::CaloID::DET_HGCAL_ENDCAP),
hitsAndFractions,
reco::CaloCluster::particleFlow, // algoID (copying from output of PFECALCSuperClusterProducer)
hitsAndFractions.at(0).first // seedId : TODO figure out if needed or if we need to put the highest energy hit
);
caloCluster.setCorrectedEnergy(emTrackster.regressed_energy());
}

edm::OrphanHandle<reco::CaloClusterCollection> caloClustersEM_h = iEvent.put(std::move(caloClustersEM));

// Fill reco::SuperCluster collection
assert(ticlSuperclusters.size() == ticlSuperclusterLinks.size());
for (std::size_t sc_i = 0; sc_i < ticlSuperclusters.size(); sc_i++) {
ticl::Trackster const& ticlSupercluster = ticlSuperclusters[sc_i];
std::vector<unsigned int> const& superclusterLink = ticlSuperclusterLinks[sc_i];

reco::CaloClusterPtrVector trackstersEMInSupercluster;
/* TODO : for now, set as SuperCluster regressed energy the sum of regressed energies of CLUE3D EM tracksters
A regression will most likely be included here */
double regressedEnergySum = 0.;
for (unsigned int tsInSc_id : superclusterLink) {
trackstersEMInSupercluster.push_back(reco::CaloClusterPtr(caloClustersEM_h, tsInSc_id));
regressedEnergySum += emTracksters[tsInSc_id].regressed_energy();
}
reco::SuperCluster& egammaSc = egammaSuperclusters->emplace_back(
ticlSupercluster.raw_energy(),
reco::SuperCluster::Point(ticlSupercluster.barycenter()),
reco::CaloClusterPtr(caloClustersEM_h,
superclusterLink[0]), // seed (first trackster in superclusterLink is the seed)
trackstersEMInSupercluster, // clusters
0., // Epreshower
0.046, // phiwidth (TODO placeholder value for now)
0.017 // etawidth (TODO placeholder value for now)
);
egammaSc.setCorrectedEnergy(regressedEnergySum);
// correctedEnergyUncertainty is left at its default value for now
}

iEvent.put(std::move(egammaSuperclusters));
}

void EGammaSuperclusterProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
edm::ParameterSetDescription desc;
desc.add<edm::InputTag>("ticlSuperClusters", edm::InputTag("ticlTracksterLinksSuperclustering"));
desc.add<edm::InputTag>("ticlTrackstersEM", edm::InputTag("ticlTrackstersCLUE3DHigh"))
->setComment("The trackster collection used before superclustering, ie CLUE3D EM tracksters");
desc.add<edm::InputTag>("layerClusters", edm::InputTag("hgcalMergeLayerClusters"))
->setComment("The layer cluster collection that goes with ticlTrackstersEM");

descriptions.add("ticlEGammaSuperClusterProducer", desc);
}

#include "FWCore/Framework/interface/MakerMacros.h"
DEFINE_FWK_MODULE(EGammaSuperclusterProducer);
8 changes: 7 additions & 1 deletion RecoHGCal/TICL/plugins/PatternRecognitionbyCLUE3D.cc
Original file line number Diff line number Diff line change
Expand Up @@ -347,7 +347,13 @@ void PatternRecognitionbyCLUE3D<TILES>::makeTracksters(
input.layerClusters,
input.layerClustersTime,
rhtools_.getPositionLayer(rhtools_.lastLayerEE(false), false).z(),
computeLocalTime_);
computeLocalTime_,
true, // energy weighting
rhtools_,
10, // min layer
10, // max layer
true // cleaning
);

if (PatternRecognitionAlgoBaseT<TILES>::algo_verbosity_ > VerbosityLevel::Advanced) {
for (auto const &t : result) {
Expand Down
Loading

0 comments on commit ac56f5d

Please sign in to comment.