6 #include <unordered_map> 10 : threshold_(conf.getParameter<double>(
"shape_threshold")),
11 distance_(conf.getParameter<double>(
"shape_distance")) {}
19 for (
const auto& energy_X : energy_X_tc) {
20 X_sum += energy_X.first * energy_X.second;
21 Etot += energy_X.first;
26 X_mean = X_sum / Etot;
31 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.
constituents();
35 for (
const auto& id_clu : clustersPtrs) {
36 if (!
pass(*id_clu.second, c3d))
39 if (layer < firstLayer)
47 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.
constituents();
48 std::unordered_map<int, float> layers_pt;
51 for (
const auto& id_cluster : clustersPtrs) {
52 if (!
pass(*id_cluster.second, c3d))
55 auto itr_insert = layers_pt.emplace(layer, 0.);
56 itr_insert.first->second += id_cluster.second->pt();
57 if (itr_insert.first->second > max_pt) {
58 max_pt = itr_insert.first->second;
66 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.
constituents();
70 for (
const auto& id_clu : clustersPtrs) {
71 if (!
pass(*id_clu.second, c3d))
74 if (layer > lastLayer)
83 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.
constituents();
85 std::vector<bool>
layers(nlayers);
86 for (
const auto& id_cluster : clustersPtrs) {
87 if (!
pass(*id_cluster.second, c3d))
89 unsigned layer = triggerGeometry.
triggerLayer(id_cluster.second->detId());
90 if (layer == 0 || layer > nlayers)
92 layers[layer - 1] =
true;
96 for (
bool layer : layers) {
101 if (length > maxlength)
110 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.
constituents();
112 std::vector<double>
layers(nlayers, 0);
113 for (
const auto& id_clu : clustersPtrs) {
114 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>&
triggerCells = id_clu.second->constituents();
116 for (
const auto& id_tc : triggerCells) {
117 if (!
pass(*id_tc.second, c3d))
119 unsigned layer = triggerGeometry.
triggerLayer(id_tc.second->detId());
120 if (layer == 0 || layer > nlayers)
122 layers[layer - 1] += id_tc.second->pt();
125 std::partial_sum(layers.begin(), layers.end(), layers.begin());
126 double pt_threshold = layers.back() *
quantile;
128 for (
double pt : layers) {
129 if (
pt > pt_threshold) {
135 double pt0 = (percentile > 0 ? layers[percentile - 1] : 0.);
136 double pt1 = layers[percentile];
137 return percentile + (pt1 - pt0 > 0. ? (pt_threshold - pt0) / (pt1 - pt0) : 0.);
141 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.
constituents();
142 std::set<double> ordered_tcs;
144 for (
const auto& id_clu : clustersPtrs) {
145 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>&
triggerCells = id_clu.second->constituents();
147 for (
const auto& id_tc : triggerCells) {
148 if (!
pass(*id_tc.second, c3d))
150 ordered_tcs.emplace(id_tc.second->pt());
151 pt_sum += id_tc.second->pt();
154 double pt_threshold = pt_sum *
quantile;
155 double partial_sum = 0.;
156 double partial_sum_prev = 0.;
158 for (
auto itr = ordered_tcs.rbegin(); itr != ordered_tcs.rend(); ++itr) {
159 partial_sum_prev = partial_sum;
162 if (partial_sum > pt_threshold) {
168 (partial_sum - partial_sum_prev > 0. ? (pt_threshold - partial_sum_prev) / (partial_sum - partial_sum_prev)
173 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.
constituents();
175 std::vector<std::pair<float, float>> tc_energy_eta;
177 for (
const auto& id_clu : clustersPtrs) {
178 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>&
triggerCells = id_clu.second->constituents();
180 for (
const auto& id_tc : triggerCells) {
181 if (!
pass(*id_tc.second, c3d))
183 tc_energy_eta.emplace_back(std::make_pair(id_tc.second->energy(), id_tc.second->eta()));
187 float SeeTot =
sigmaXX(tc_energy_eta, c3d.
eta());
193 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.
constituents();
195 std::vector<std::pair<float, float>> tc_energy_phi;
197 for (
const auto& id_clu : clustersPtrs) {
198 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>&
triggerCells = id_clu.second->constituents();
200 for (
const auto& id_tc : triggerCells) {
201 if (!
pass(*id_tc.second, c3d))
203 tc_energy_phi.emplace_back(std::make_pair(id_tc.second->energy(), id_tc.second->phi()));
213 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.
constituents();
215 std::vector<std::pair<float, float>> tc_energy_r;
217 for (
const auto& id_clu : clustersPtrs) {
218 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>&
triggerCells = id_clu.second->constituents();
220 for (
const auto& id_tc : triggerCells) {
221 if (!
pass(*id_tc.second, c3d))
223 float r = (id_tc.second->position().z() != 0.
224 ?
std::sqrt(
pow(id_tc.second->position().x(), 2) +
pow(id_tc.second->position().y(), 2)) /
225 std::abs(id_tc.second->position().z())
227 tc_energy_r.emplace_back(std::make_pair(id_tc.second->energy(),
r));
231 float r_mean =
meanX(tc_energy_r);
232 float Szz =
sigmaXX(tc_energy_r, r_mean);
238 std::unordered_map<int, std::vector<std::pair<float, float>>> tc_layer_energy_eta;
239 std::unordered_map<int, LorentzVector> layer_LV;
241 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.
constituents();
243 for (
const auto& id_clu : clustersPtrs) {
246 layer_LV[layer] += id_clu.second->p4();
248 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>&
triggerCells = id_clu.second->constituents();
250 for (
const auto& id_tc : triggerCells) {
251 if (!
pass(*id_tc.second, c3d))
253 tc_layer_energy_eta[layer].emplace_back(std::make_pair(id_tc.second->energy(), id_tc.second->eta()));
257 float SigmaEtaEtaMax = 0;
259 for (
auto& tc_iter : tc_layer_energy_eta) {
260 const std::vector<std::pair<float, float>>& energy_eta_layer = tc_iter.second;
262 float SigmaEtaEtaLayer =
sigmaXX(energy_eta_layer, LV_layer.eta());
263 if (SigmaEtaEtaLayer > SigmaEtaEtaMax)
264 SigmaEtaEtaMax = SigmaEtaEtaLayer;
267 return SigmaEtaEtaMax;
271 std::unordered_map<int, std::vector<std::pair<float, float>>> tc_layer_energy_phi;
272 std::unordered_map<int, LorentzVector> layer_LV;
274 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.
constituents();
276 for (
const auto& id_clu : clustersPtrs) {
279 layer_LV[layer] += id_clu.second->p4();
281 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>&
triggerCells = id_clu.second->constituents();
283 for (
const auto& id_tc : triggerCells) {
284 if (!
pass(*id_tc.second, c3d))
286 tc_layer_energy_phi[layer].emplace_back(std::make_pair(id_tc.second->energy(), id_tc.second->phi()));
290 float SigmaPhiPhiMax = 0;
292 for (
auto& tc_iter : tc_layer_energy_phi) {
293 const std::vector<std::pair<float, float>>& energy_phi_layer = tc_iter.second;
295 float SigmaPhiPhiLayer =
sigmaPhiPhi(energy_phi_layer, LV_layer.phi());
296 if (SigmaPhiPhiLayer > SigmaPhiPhiMax)
297 SigmaPhiPhiMax = SigmaPhiPhiLayer;
300 return SigmaPhiPhiMax;
304 std::unordered_map<int, std::vector<std::pair<float, float>>> tc_layer_energy_r;
306 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.
constituents();
308 for (
const auto& id_clu : clustersPtrs) {
311 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>&
triggerCells = id_clu.second->constituents();
313 for (
const auto& id_tc : triggerCells) {
314 if (!
pass(*id_tc.second, c3d))
316 float r = (id_tc.second->position().z() != 0.
317 ?
std::sqrt(
pow(id_tc.second->position().x(), 2) +
pow(id_tc.second->position().y(), 2)) /
318 std::abs(id_tc.second->position().z())
320 tc_layer_energy_r[layer].emplace_back(std::make_pair(id_tc.second->energy(),
r));
324 float SigmaRRMax = 0;
326 for (
auto& tc_iter : tc_layer_energy_r) {
327 const std::vector<std::pair<float, float>>& energy_r_layer = tc_iter.second;
328 float r_mean_layer =
meanX(energy_r_layer);
329 float SigmaRRLayer =
sigmaXX(energy_r_layer, r_mean_layer);
330 if (SigmaRRLayer > SigmaRRMax)
331 SigmaRRMax = SigmaRRLayer;
338 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.
constituents();
340 std::unordered_map<int, std::vector<edm::Ptr<l1t::HGCalTriggerCell>>> layers_tcs;
341 for (
const auto& id_clu : clustersPtrs) {
343 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>&
triggerCells = id_clu.second->constituents();
344 for (
const auto& id_tc : triggerCells) {
345 if (!
pass(*id_tc.second, c3d))
347 layers_tcs[layer].emplace_back(id_tc.second);
352 std::unordered_map<int, std::vector<std::pair<float, float>>> tc_layers_energy_r;
353 for (
const auto& layer_tcs : layers_tcs) {
354 int layer = layer_tcs.first;
356 for (
const auto& tc : layer_tcs.second) {
357 if (tc->energy() > max_tc->
energy())
360 for (
const auto& tc : layer_tcs.second) {
361 double dx = tc->position().x() - max_tc->
position().
x();
362 double dy = tc->position().y() - max_tc->
position().
y();
363 double distance_to_max =
std::sqrt(dx * dx + dy * dy);
364 if (distance_to_max < radius) {
365 float r = (tc->position().z() != 0.
366 ?
std::sqrt(tc->position().x() * tc->position().x() + tc->position().y() * tc->position().y()) /
369 tc_layers_energy_r[layer].emplace_back(std::make_pair(tc->energy(),
r));
375 std::vector<std::pair<float, float>> layers_energy_srr2;
376 for (
const auto& layer_energy_r : tc_layers_energy_r) {
377 const auto& energies_r = layer_energy_r.second;
378 float r_mean_layer =
meanX(energies_r);
379 float srr =
sigmaXX(energies_r, r_mean_layer);
380 double energy_sum = 0.;
381 for (
const auto& energy_r : energies_r) {
382 energy_sum += energy_r.first;
384 layers_energy_srr2.emplace_back(std::make_pair(energy_sum, srr * srr));
387 float srr2_mean =
meanX(layers_energy_srr2);
392 std::unordered_map<int, float> layer_energy;
394 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.
constituents();
396 for (
const auto& id_clu : clustersPtrs) {
397 if (!
pass(*id_clu.second, c3d))
400 layer_energy[layer] += id_clu.second->energy();
405 for (
const auto& layer : layer_energy) {
406 if (layer.second > EMax)
414 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.
constituents();
416 std::vector<std::pair<float, float>> tc_energy_z;
418 for (
const auto& id_clu : clustersPtrs) {
419 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>&
triggerCells = id_clu.second->constituents();
421 for (
const auto& id_tc : triggerCells) {
422 if (!
pass(*id_tc.second, c3d))
424 tc_energy_z.emplace_back(id_tc.second->energy(), id_tc.second->position().z());
428 return meanX(tc_energy_z);
432 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.
constituents();
434 std::vector<std::pair<float, float>> tc_energy_z;
436 for (
const auto& id_clu : clustersPtrs) {
437 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>&
triggerCells = id_clu.second->constituents();
439 for (
const auto& id_tc : triggerCells) {
440 if (!
pass(*id_tc.second, c3d))
442 tc_energy_z.emplace_back(std::make_pair(id_tc.second->energy(), id_tc.second->position().z()));
446 float z_mean =
meanX(tc_energy_z);
447 float Szz =
sigmaXX(tc_energy_z, z_mean);
453 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& cellsPtrs = c2d.
constituents();
455 std::vector<std::pair<float, float>> tc_energy_eta;
457 for (
const auto& id_cell : cellsPtrs) {
458 if (!
pass(*id_cell.second, c2d))
460 tc_energy_eta.emplace_back(std::make_pair(id_cell.second->energy(), id_cell.second->eta()));
463 float See =
sigmaXX(tc_energy_eta, c2d.
eta());
469 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& cellsPtrs = c2d.
constituents();
471 std::vector<std::pair<float, float>> tc_energy_phi;
473 for (
const auto& id_cell : cellsPtrs) {
474 if (!
pass(*id_cell.second, c2d))
476 tc_energy_phi.emplace_back(std::make_pair(id_cell.second->energy(), id_cell.second->phi()));
485 const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& cellsPtrs = c2d.
constituents();
487 std::vector<std::pair<float, float>> tc_energy_r;
489 for (
const auto& id_cell : cellsPtrs) {
490 if (!
pass(*id_cell.second, c2d))
492 float r = (id_cell.second->position().z() != 0.
493 ?
std::sqrt(
pow(id_cell.second->position().x(), 2) +
pow(id_cell.second->position().y(), 2)) /
494 std::abs(id_cell.second->position().z())
496 tc_energy_r.emplace_back(std::make_pair(id_cell.second->energy(),
r));
499 float r_mean =
meanX(tc_energy_r);
500 float Srr =
sigmaXX(tc_energy_r, r_mean);
float meanZ(const l1t::HGCalMulticluster &c3d) const
float layer10percent() const
float sigmaRRMax(const l1t::HGCalMulticluster &c3d) const
double eta() const final
momentum pseudorapidity
bool pass(const T &obj, const Tref &ref) const
float sigmaPhiPhiTot() const
float triggerCells90percent() const
float sigmaXX(const std::vector< pair< float, float >> &energy_X_tc, const float X_cluster) const
float sigmaEtaEtaMax() const
float sigmaPhiPhi(const std::vector< pair< float, float >> &energy_phi_tc, const float phi_cluster) const
int coreShowerLength(const l1t::HGCalMulticluster &c3d, const HGCalTriggerGeometryBase &triggerGeometry) const
float sigmaRRTot(const l1t::HGCalMulticluster &c3d) const
HGCalTriggerTools triggerTools_
float sigmaPhiPhiTot(const l1t::HGCalMulticluster &c3d) const
int firstLayer(const l1t::HGCalMulticluster &c3d) const
int lastLayer(const l1t::HGCalMulticluster &c3d) const
float percentileLayer(const l1t::HGCalMulticluster &c3d, const HGCalTriggerGeometryBase &triggerGeometry, float quantile=0.5) const
const GlobalPoint & position() const
float sigmaRRMean(const l1t::HGCalMulticluster &c3d, float radius=5.) const
float layer90percent() const
double energy() const final
energy
float sigmaEtaEtaMax(const l1t::HGCalMulticluster &c3d) const
float percentileTriggerCells(const l1t::HGCalMulticluster &c3d, float quantile=0.5) const
Abs< T >::type abs(const T &t)
virtual unsigned triggerLayer(const unsigned id) const =0
math::XYZTLorentzVector LorentzVector
float meanX(const std::vector< pair< float, float >> &energy_X_tc) const
float sigmaEtaEtaTot() const
int coreShowerLength() const
float layer50percent() const
void fillShapes(l1t::HGCalMulticluster &, const HGCalTriggerGeometryBase &) const
int showerLength(const l1t::HGCalMulticluster &c3d) const
float sigmaRRMean() const
float triggerCells67percent() const
float sigmaZZ(const l1t::HGCalMulticluster &c3d) const
float sigmaPhiPhiMax(const l1t::HGCalMulticluster &c3d) const
int maxLayer(const l1t::HGCalMulticluster &c3d) const
float zBarycenter() const
float sigmaPhiPhiMax() const
float eMax(const l1t::HGCalMulticluster &c3d) const
double phi() const final
momentum azimuthal angle
float sigmaEtaEtaTot(const l1t::HGCalMulticluster &c3d) const
Power< A, B >::type pow(const A &a, const B &b)
const std::unordered_map< uint32_t, edm::Ptr< C > > & constituents() const