LDMX Software
ecal::EcalClusterProducer Class Reference

Simple algorithm that does clustering in the ECal. More...

#include <EcalClusterProducer.h>

Public Member Functions

 EcalClusterProducer (const std::string &name, framework::Process &process)
 
void configure (framework::config::Parameters &parameters) override
 Configure the processor using the given user specified parameters.
 
void produce (framework::Event &event) override
 Process the event and put new data products into it.
 
- Public Member Functions inherited from framework::Producer
 Producer (const std::string &name, Process &process)
 Class constructor.
 
virtual void process (Event &event) final
 Processing an event for a Producer is calling produce.
 
- Public Member Functions inherited from framework::EventProcessor
 DECLARE_FACTORY (EventProcessor, EventProcessor *, const std::string &, Process &)
 declare that we have a factory for this class
 
 EventProcessor (const std::string &name, Process &process)
 Class constructor.
 
virtual ~EventProcessor ()=default
 Class destructor.
 
virtual void beforeNewRun (ldmx::RunHeader &run_header)
 Callback for Producers to add parameters to the run header before conditions are initialized.
 
virtual void onNewRun (const ldmx::RunHeader &run_header)
 Callback for the EventProcessor to take any necessary action when the run being processed changes.
 
virtual void onFileOpen (EventFile &event_file)
 Callback for the EventProcessor to take any necessary action when a new event input ROOT file is opened.
 
virtual void onFileClose (EventFile &event_file)
 Callback for the EventProcessor to take any necessary action when a event input ROOT file is closed.
 
virtual void onProcessStart ()
 Callback for the EventProcessor to take any necessary action when the processing of events starts, such as creating histograms.
 
virtual void onProcessEnd ()
 Callback for the EventProcessor to take any necessary action when the processing of events finishes, such as calculating job-summary quantities.
 
template<class T >
const T & getCondition (const std::string &condition_name)
 Access a conditions object for the current event.
 
TDirectory * getHistoDirectory ()
 Access/create a directory in the histogram file for this event processor to create histograms and analysis tuples.
 
void setStorageHint (framework::StorageControl::Hint hint)
 Mark the current event as having the given storage control hint from this module_.
 
void setStorageHint (framework::StorageControl::Hint hint, const std::string &purposeString)
 Mark the current event as having the given storage control hint from this module and the given purpose string.
 
int getLogFrequency () const
 Get the current logging frequency from the process.
 
int getRunNumber () const
 Get the run number from the process.
 
std::string getName () const
 Get the processor name.
 
void createHistograms (const std::vector< framework::config::Parameters > &histos)
 Internal function which is used to create histograms passed from the python configuration @parma histos vector of Parameters that configure histograms to create.
 

Private Attributes

double seed_threshold_ {0}
 
double cutoff_ {0}
 
double dc_ {0}
 
double rhoc_ {0}
 
double deltac_ {0}
 
double deltao_ {0}
 
float min_hit_energy_ {1.}
 
std::string rec_hit_coll_name_
 
std::string rec_hit_pass_name_
 
std::string algo_coll_name_
 
std::string cluster_coll_name_
 
bool clue_
 
int nbr_of_layers_
 
bool reclustering_
 
std::string algorithm_
 
CLUE3D::Params clue3d_
 
TString algo_name_
 The name of the cluster algorithm used.
 

Additional Inherited Members

- Protected Member Functions inherited from framework::EventProcessor
void abortEvent ()
 Abort the event immediately.
 
- Protected Attributes inherited from framework::EventProcessor
HistogramPool histograms_
 helper object for making and filling histograms
 
NtupleManager & ntuple_ {NtupleManager::getInstance()}
 Manager for any ntuples.
 
logging::logger the_log_
 The logger for this EventProcessor.
 

Detailed Description

Simple algorithm that does clustering in the ECal.

Definition at line 23 of file EcalClusterProducer.h.

Constructor & Destructor Documentation

◆ EcalClusterProducer()

ecal::EcalClusterProducer::EcalClusterProducer ( const std::string & name,
framework::Process & process )

Definition at line 20 of file EcalClusterProducer.cxx.

22 : Producer(name, process) {}
Producer(const std::string &name, Process &process)
Class constructor.
virtual void process(Event &event) final
Processing an event for a Producer is calling produce.

Member Function Documentation

◆ configure()

void ecal::EcalClusterProducer::configure ( framework::config::Parameters & parameters)
overridevirtual

Configure the processor using the given user specified parameters.

Parameters
parametersSet of parameters used to configure this processor.

Reimplemented from framework::EventProcessor.

Definition at line 24 of file EcalClusterProducer.cxx.

24 {
25 cutoff_ = ps.get<double>("cutoff");
26 seed_threshold_ = ps.get<double>("seed_threshold");
27
28 dc_ = ps.get<double>("dc");
29 rhoc_ = ps.get<double>("rhoc");
30 deltac_ = ps.get<double>("deltac");
31 deltao_ = ps.get<double>("deltao");
32
33 rec_hit_coll_name_ = ps.get<std::string>("rec_hit_coll_name");
34 rec_hit_pass_name_ = ps.get<std::string>("rec_hit_pass_name");
35 algo_coll_name_ = ps.get<std::string>("algo_coll_name");
36 algo_name_ = ps.get<std::string>("algo_name");
37 cluster_coll_name_ = ps.get<std::string>("cluster_coll_name");
38 clue_ = ps.get<bool>("clue");
39 nbr_of_layers_ = ps.get<int>("nbr_of_layers");
40 reclustering_ = ps.get<bool>("reclustering");
41 algorithm_ = ps.get<std::string>("algorithm", "clue");
42
43 clue3d_.hit_ecut_ = ps.get<double>("clue3d_hit_ecut", clue3d_.hit_ecut_);
44 clue3d_.dc2d_ = ps.get<double>("clue3d_dc2d", clue3d_.dc2d_);
45 clue3d_.kernel2d_ = ps.get<double>("clue3d_kernel2d", clue3d_.kernel2d_);
46 clue3d_.deltas2d_ = ps.get<double>("clue3d_deltas2d", clue3d_.deltas2d_);
47 clue3d_.deltao2d_ = ps.get<double>("clue3d_deltao2d", clue3d_.deltao2d_);
48 clue3d_.rhoc2d_ = ps.get<double>("clue3d_rhoc2d", clue3d_.rhoc2d_);
49 clue3d_.sibling_layers_ =
50 ps.get<int>("clue3d_sibling_layers", clue3d_.sibling_layers_);
51 clue3d_.dens_xy3d_ = ps.get<double>("clue3d_dens_xy3d", clue3d_.dens_xy3d_);
52 clue3d_.kernel3d_ = ps.get<double>("clue3d_kernel3d", clue3d_.kernel3d_);
53 clue3d_.crit_xy3d_ = ps.get<double>("clue3d_crit_xy3d", clue3d_.crit_xy3d_);
54 clue3d_.crit_rho3d_ =
55 ps.get<double>("clue3d_crit_rho3d", clue3d_.crit_rho3d_);
56 clue3d_.self_frac3d_ =
57 ps.get<double>("clue3d_self_frac3d", clue3d_.self_frac3d_);
58 clue3d_.outlier_mult3d_ =
59 ps.get<double>("clue3d_outlier_mult3d", clue3d_.outlier_mult3d_);
60 clue3d_.attach_radius_ =
61 ps.get<double>("clue3d_attach_radius", clue3d_.attach_radius_);
62 clue3d_.sc_ratio_ = ps.get<double>("clue3d_sc_ratio", clue3d_.sc_ratio_);
63 clue3d_.sc_radius_ = ps.get<double>("clue3d_sc_radius", clue3d_.sc_radius_);
64 clue3d_.sc_min_layers_ =
65 ps.get<int>("clue3d_sc_min_layers", clue3d_.sc_min_layers_);
66 clue3d_.coincide_dist_ =
67 ps.get<double>("clue3d_coincide_dist", clue3d_.coincide_dist_);
68 clue3d_.ghost_ratio_ =
69 ps.get<double>("clue3d_ghost_ratio", clue3d_.ghost_ratio_);
70 clue3d_.ghost_dist_ =
71 ps.get<double>("clue3d_ghost_dist", clue3d_.ghost_dist_);
72 clue3d_.min_cluster_energy_ =
73 ps.get<double>("clue3d_min_cluster_energy", clue3d_.min_cluster_energy_);
74 clue3d_.share_ = ps.get<bool>("clue3d_share", clue3d_.share_);
75 clue3d_.share_iters_ =
76 ps.get<int>("clue3d_share_iters", clue3d_.share_iters_);
77 clue3d_.share_min_energy_ =
78 ps.get<double>("clue3d_share_min_energy", clue3d_.share_min_energy_);
79 clue3d_.profile_lambda_ = ps.get<std::vector<double>>(
80 "clue3d_profile_lambda", clue3d_.profile_lambda_);
81 clue3d_.w0_ = ps.get<double>("clue3d_w0", clue3d_.w0_);
82 clue3d_.pos_window_ =
83 ps.get<double>("clue3d_pos_window", clue3d_.pos_window_);
84}
TString algo_name_
The name of the cluster algorithm used.

References framework::config::Parameters::get().

◆ produce()

void ecal::EcalClusterProducer::produce ( framework::Event & event)
overridevirtual

Process the event and put new data products into it.

Parameters
eventThe Event to process.

Implements framework::Producer.

Definition at line 86 of file EcalClusterProducer.cxx.

86 {
87 const auto& ecal_hits{event.getCollection<ldmx::EcalHit>(rec_hit_coll_name_,
88 rec_hit_pass_name_)};
89 if (ecal_hits.size() == 0) {
90 // keep the collection present for downstream readers
91 ldmx_log(debug) << "No ECal hits found";
92 std::vector<ldmx::EcalCluster> empty;
93 event.add(cluster_coll_name_, empty);
94 return;
95 }
96
97 if (algorithm_ == "clue3d") {
98 CLUE3D algo(clue3d_);
99 auto clusters = algo.cluster(ecal_hits);
100 std::vector<ldmx::EcalCluster> ecal_clusters;
101 for (const auto& c : clusters) {
102 ldmx::EcalCluster cluster;
103 cluster.setEnergy(c.e_);
104 cluster.setCentroidXYZ(c.x_, c.y_, c.z_);
105 cluster.setDXDZ(c.dxdz_);
106 cluster.setDYDZ(c.dydz_);
107 cluster.setRMSXYZ(c.rms_x_, c.rms_y_, c.rms_z_);
108 cluster.setFirstLayerCentroidXYZ(c.x0_, c.y0_, c.z0_);
109 // 3D clusters span all layers
110 cluster.setLayer(0);
111 std::vector<unsigned int> ids;
112 std::vector<float> hx, hy, hz, he;
113 for (std::size_t i = 0; i < c.hits_.size(); i++) {
114 const auto& h = ecal_hits[c.hits_[i]];
115 ids.push_back(h.getID());
116 hx.push_back(h.getXPos());
117 hy.push_back(h.getYPos());
118 hz.push_back(h.getZPos());
119 // shared energy of the hit
120 he.push_back(c.frac_[i] * h.getEnergy());
121 }
122 cluster.setNHits(ids.size());
123 cluster.setIDs(ids);
124 cluster.setHitValsX(hx);
125 cluster.setHitValsY(hy);
126 cluster.setHitValsZ(hz);
127 cluster.setHitValsE(he);
128 ecal_clusters.push_back(cluster);
129 }
130 event.add(cluster_coll_name_, ecal_clusters);
131 } else if (clue_) {
132 ldmx_log(info) << "Using CLUE clustering algorithm";
133 CLUE cf;
134 cf.cluster(ecal_hits, dc_, rhoc_, deltac_, deltao_, nbr_of_layers_,
135 reclustering_);
136 ldmx_log(debug) << "CLUE algorithm finished clustering";
137 std::vector<IntermediateCluster> interm_cluster = cf.getClusters();
138 std::vector<IntermediateCluster> f_interm_cluster =
139 cf.getFirstLayerCentroids();
140 ldmx_log(debug) << "Got " << interm_cluster.size() << " clusters with "
141 << f_interm_cluster.size() << " first layer centroids";
142
143 auto n_loops = cf.getNLoops();
144 ldmx_log(debug) << "Number of clustering loops: " << n_loops;
145
146 if (reclustering_) {
147 ldmx_log(debug) << "Reclustererd initial number of clusters: "
148 << cf.getInitialClusterNbr()
149 << ", final number of clusters: "
150 << interm_cluster.size();
151 }
152
153 std::vector<ldmx::EcalCluster> ecal_clusters;
154 ldmx_log(debug) << "Filling " << interm_cluster.size()
155 << " clusters into ecal_clusters";
156 for (size_t cluster_indx = 0; cluster_indx < interm_cluster.size();
157 cluster_indx++) {
158 ldmx::EcalCluster cluster;
159
160 cluster.setEnergy(interm_cluster[cluster_indx].energy());
161 cluster.setCentroidXYZ(interm_cluster[cluster_indx].centroidX(),
162 interm_cluster[cluster_indx].centroidY(),
163 interm_cluster[cluster_indx].centroidZ());
164 cluster.setFirstLayerCentroidXYZ(
165 f_interm_cluster[cluster_indx].centroidX(),
166 f_interm_cluster[cluster_indx].centroidY(),
167 f_interm_cluster[cluster_indx].centroidZ());
168 cluster.setNHits(interm_cluster[cluster_indx].hits().size());
169 cluster.addHits(interm_cluster[cluster_indx].hits());
170 cluster.addFirstLayerHits(f_interm_cluster[cluster_indx].hits());
171
172 float cl_x(0), cl_y(0), cl_z(0), cl_xx(0), cl_yy(0), cl_zz(0);
173 float cl_w = 1; // weight
174 float sumw = 0;
175
176 for (auto hit : interm_cluster[cluster_indx].hits()) {
177 if (hit->getEnergy() < min_hit_energy_) continue;
178 cl_w = log(hit->getEnergy()) - log(min_hit_energy_);
179 cl_x += cl_w * hit->getXPos();
180 cl_y += cl_w * hit->getYPos();
181 cl_z += cl_w * hit->getZPos();
182 cl_xx += cl_w * hit->getXPos() * hit->getXPos();
183 cl_yy += cl_w * hit->getYPos() * hit->getYPos();
184 cl_zz += cl_w * hit->getZPos() * hit->getZPos();
185 sumw += cl_w;
186 } // over hits
187 // could probably get this as cluster.getCentroidX() instead
188 cl_x /= sumw; // now is <x_>
189 cl_y /= sumw;
190 cl_z /= sumw;
191 cl_xx /= sumw; // now is <x_^2>
192 cl_yy /= sumw;
193 cl_zz /= sumw;
194 cl_xx = sqrt(cl_xx - cl_x * cl_x); // now is sqrt(<x_^2>-<x_>^2)
195 cl_yy = sqrt(cl_yy - cl_y * cl_y);
196 cl_zz = sqrt(cl_zz - cl_z * cl_z);
197
198 cluster.setRMSXYZ(cl_xx, cl_yy, cl_zz);
199 cluster.setLayer(interm_cluster[cluster_indx].layer());
200 ldmx_log(trace) << "Cluster " << cluster_indx
201 << " energy: " << cluster.getEnergy()
202 << ", nHits: " << cluster.getNHits() << ", centroid: ("
203 << cluster.getCentroidX() << ", "
204 << cluster.getCentroidY() << ", "
205 << cluster.getCentroidZ() << "), " << "RMS : ("
206 << cluster.getRMSX() << ", " << cluster.getRMSY() << ", "
207 << cluster.getRMSZ() << ")" << ", first layer centroid: ("
208 << cluster.getFirstLayerCentroidX() << ", "
209 << cluster.getFirstLayerCentroidY() << ", "
210 << cluster.getFirstLayerCentroidZ() << ")";
211
212 ecal_clusters.push_back(cluster);
213 }
214 ldmx_log(debug) << "Filled " << ecal_clusters.size()
215 << " clusters into ecal_clusters";
216 event.add(cluster_coll_name_, ecal_clusters);
217 } else {
218 ldmx_log(info) << "Using simple clustering algorithm " << algo_name_;
220
221 for (const ldmx::EcalHit& hit : ecal_hits) {
222 // Skip zero energy digis.
223 if (hit.getEnergy() == 0) {
224 continue;
225 }
226 cf.add(hit);
227 }
228
229 cf.cluster(seed_threshold_, cutoff_);
230 auto interm_cluster = cf.getClusters();
231 std::map<int, double> c_weights = cf.getWeights();
232
233 ldmx::ClusterAlgoResult algo_result;
234 algo_result.set(algo_name_, 3, c_weights.rbegin()->first);
235 algo_result.setAlgoVar(0, cutoff_);
236 algo_result.setAlgoVar(1, seed_threshold_);
237 algo_result.setAlgoVar(2, cf.getNSeeds());
238
239 std::map<int, double>::iterator it = c_weights.begin();
240 for (it = c_weights.begin(); it != c_weights.end(); it++) {
241 algo_result.setWeight(it->first, it->second / 100);
242 }
243
244 std::vector<ldmx::EcalCluster> ecal_clusters;
245 for (size_t cluster_indx = 0; cluster_indx < interm_cluster.size();
246 cluster_indx++) {
247 ldmx::EcalCluster cluster;
248
249 cluster.setEnergy(interm_cluster[cluster_indx].energy());
250 cluster.setCentroidXYZ(interm_cluster[cluster_indx].centroidX(),
251 interm_cluster[cluster_indx].centroidY(),
252 interm_cluster[cluster_indx].centroidZ());
253 cluster.setNHits(interm_cluster[cluster_indx].hits().size());
254 cluster.addHits(interm_cluster[cluster_indx].hits());
255 ecal_clusters.push_back(cluster);
256 }
257
258 event.add(cluster_coll_name_, ecal_clusters);
259 event.add(algo_coll_name_, algo_result);
260 } // end on simple clustering
261}
double getRMSX() const
rms in x
void setNHits(int nHits)
Sets total number of hits in the cluster.
Definition CaloCluster.h:64
double getCentroidZ() const
centroid z-location
double getCentroidX() const
centroid x-location
void setCentroidXYZ(double centroid_x, double centroid_y, double centroid_z)
Sets the three coordinates of the cluster centroid.
Definition CaloCluster.h:84
double getRMSZ() const
rms in z
double getRMSY() const
rms in y
void setEnergy(double energy)
Sets total energy for the cluster.
Definition CaloCluster.h:58
void setIDs(std::vector< unsigned int > &hitIDs)
Sets a sorted vector for the IDs of the hits_ that make up the cluster.
Definition CaloCluster.h:71
double getCentroidY() const
centroid y-location
Contains details about the clustering algorithm.
void setAlgoVar(int element, double value)
Set an algorithm variable.
void set(const TString &name, int nvar)
Set name and number of variables of cluster algo.
void setWeight(int nClusters, double weight)
Set a weight when number of clusters reached.
Stores cluster information from the ECal.
Definition EcalCluster.h:20
void addHits(const std::vector< const ldmx::EcalHit * > &hits)
Take in the hits that make up the cluster.
Stores reconstructed hit information from the ECAL.
Definition EcalHit.h:19
A templated agglomerative clustering algorithm.
void cluster(double seed_threshold, double cutoff)
Run the clustering algorithm.
std::vector< ClusterType > getClusters() const
Get the final clusters after filtering by seed threshold.
std::map< int, double > getWeights() const
Get the transition weights (cluster count -> minimum weight at that step).
int getNSeeds() const
Get the number of seed clusters found.
void add(const HitType &hit)
Add a hit to be clustered using its stored position.

References recon::TemplatedClusterFinder< HitType, WeightClass >::add(), ldmx::EcalCluster::addHits(), ecal::CLUE3D::cluster(), recon::TemplatedClusterFinder< HitType, WeightClass >::cluster(), ldmx::CaloCluster::getCentroidX(), ldmx::CaloCluster::getCentroidY(), ldmx::CaloCluster::getCentroidZ(), recon::TemplatedClusterFinder< HitType, WeightClass >::getClusters(), recon::TemplatedClusterFinder< HitType, WeightClass >::getNSeeds(), ldmx::CaloCluster::getRMSX(), ldmx::CaloCluster::getRMSY(), ldmx::CaloCluster::getRMSZ(), recon::TemplatedClusterFinder< HitType, WeightClass >::getWeights(), ldmx::ClusterAlgoResult::set(), ldmx::ClusterAlgoResult::setAlgoVar(), ldmx::CaloCluster::setCentroidXYZ(), ldmx::CaloCluster::setEnergy(), ldmx::CaloCluster::setIDs(), ldmx::CaloCluster::setNHits(), and ldmx::ClusterAlgoResult::setWeight().

Member Data Documentation

◆ algo_coll_name_

std::string ecal::EcalClusterProducer::algo_coll_name_
private

Definition at line 50 of file EcalClusterProducer.h.

◆ algo_name_

TString ecal::EcalClusterProducer::algo_name_
private

The name of the cluster algorithm used.

Definition at line 62 of file EcalClusterProducer.h.

◆ algorithm_

std::string ecal::EcalClusterProducer::algorithm_
private

Definition at line 58 of file EcalClusterProducer.h.

◆ clue3d_

CLUE3D::Params ecal::EcalClusterProducer::clue3d_
private

Definition at line 59 of file EcalClusterProducer.h.

◆ clue_

bool ecal::EcalClusterProducer::clue_
private

Definition at line 53 of file EcalClusterProducer.h.

◆ cluster_coll_name_

std::string ecal::EcalClusterProducer::cluster_coll_name_
private

Definition at line 51 of file EcalClusterProducer.h.

◆ cutoff_

double ecal::EcalClusterProducer::cutoff_ {0}
private

Definition at line 39 of file EcalClusterProducer.h.

39{0};

◆ dc_

double ecal::EcalClusterProducer::dc_ {0}
private

Definition at line 41 of file EcalClusterProducer.h.

41{0};

◆ deltac_

double ecal::EcalClusterProducer::deltac_ {0}
private

Definition at line 43 of file EcalClusterProducer.h.

43{0};

◆ deltao_

double ecal::EcalClusterProducer::deltao_ {0}
private

Definition at line 44 of file EcalClusterProducer.h.

44{0};

◆ min_hit_energy_

float ecal::EcalClusterProducer::min_hit_energy_ {1.}
private

Definition at line 46 of file EcalClusterProducer.h.

46{1.};

◆ nbr_of_layers_

int ecal::EcalClusterProducer::nbr_of_layers_
private

Definition at line 54 of file EcalClusterProducer.h.

◆ rec_hit_coll_name_

std::string ecal::EcalClusterProducer::rec_hit_coll_name_
private

Definition at line 48 of file EcalClusterProducer.h.

◆ rec_hit_pass_name_

std::string ecal::EcalClusterProducer::rec_hit_pass_name_
private

Definition at line 49 of file EcalClusterProducer.h.

◆ reclustering_

bool ecal::EcalClusterProducer::reclustering_
private

Definition at line 55 of file EcalClusterProducer.h.

◆ rhoc_

double ecal::EcalClusterProducer::rhoc_ {0}
private

Definition at line 42 of file EcalClusterProducer.h.

42{0};

◆ seed_threshold_

double ecal::EcalClusterProducer::seed_threshold_ {0}
private

Definition at line 38 of file EcalClusterProducer.h.

38{0};

The documentation for this class was generated from the following files: