src/backend/HGCalShowerShape.cc

0001 #include "L1Trigger/L1THGCal/interface/backend/HGCalShowerShape.h"
0002 #include "DataFormats/L1THGCal/interface/HGCalMulticluster.h"
0003 #include "DataFormats/L1THGCal/interface/HGCalCluster.h"
0004 #include "DataFormats/Math/interface/deltaPhi.h"
0005
0006 #include <unordered_map>
0007 #include <numeric>
0008
0009 HGCalShowerShape::HGCalShowerShape(const edm::ParameterSet& conf)
0010     : threshold_(conf.getParameter<double>("shape_threshold")),
0011       distance_(conf.getParameter<double>("shape_distance")) {}
0012
0013 //Compute energy-weighted mean of any variable X in the cluster
0014
0015 float HGCalShowerShape::meanX(const std::vector<pair<float, float>>& energy_X_tc) const {
0016   float Etot = 0;
0017   float X_sum = 0;
0018
0019   for (const auto& energy_X : energy_X_tc) {
0020     X_sum += energy_X.first * energy_X.second;
0021     Etot += energy_X.first;
0022   }
0023
0024   float X_mean = 0;
0025   if (Etot > 0)
0026     X_mean = X_sum / Etot;
0027   return X_mean;
0028 }
0029
0030 int HGCalShowerShape::firstLayer(const l1t::HGCalMulticluster& c3d) const {
0031   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0032
0033   int firstLayer = 999;
0034
0035   for (const auto& id_clu : clustersPtrs) {
0036     if (!pass(*id_clu.second, c3d))
0037       continue;
0038     int layer = triggerTools_.layerWithOffset(id_clu.second->detId());
0039     if (layer < firstLayer)
0040       firstLayer = layer;
0041   }
0042
0043   return firstLayer;
0044 }
0045
0046 int HGCalShowerShape::maxLayer(const l1t::HGCalMulticluster& c3d) const {
0047   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0048   std::unordered_map<int, float> layers_pt;
0049   float max_pt = 0.;
0050   int max_layer = 0;
0051   for (const auto& id_cluster : clustersPtrs) {
0052     if (!pass(*id_cluster.second, c3d))
0053       continue;
0054     unsigned layer = triggerTools_.layerWithOffset(id_cluster.second->detId());
0055     auto itr_insert = layers_pt.emplace(layer, 0.);
0056     itr_insert.first->second += id_cluster.second->pt();
0057     if (itr_insert.first->second > max_pt) {
0058       max_pt = itr_insert.first->second;
0059       max_layer = layer;
0060     }
0061   }
0062   return max_layer;
0063 }
0064
0065 int HGCalShowerShape::lastLayer(const l1t::HGCalMulticluster& c3d) const {
0066   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0067
0068   int lastLayer = -999;
0069
0070   for (const auto& id_clu : clustersPtrs) {
0071     if (!pass(*id_clu.second, c3d))
0072       continue;
0073     int layer = triggerTools_.layerWithOffset(id_clu.second->detId());
0074     if (layer > lastLayer)
0075       lastLayer = layer;
0076   }
0077
0078   return lastLayer;
0079 }
0080
0081 int HGCalShowerShape::coreShowerLength(const l1t::HGCalMulticluster& c3d,
0082                                        const HGCalTriggerGeometryBase& triggerGeometry) const {
0083   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0084   unsigned nlayers = triggerTools_.layers(ForwardSubdetector::ForwardEmpty);
0085   std::vector<bool> layers(nlayers);
0086   for (const auto& id_cluster : clustersPtrs) {
0087     if (!pass(*id_cluster.second, c3d))
0088       continue;
0089     unsigned layer = triggerGeometry.triggerLayer(id_cluster.second->detId());
0090     if (triggerTools_.isNose(id_cluster.second->detId()))
0091       nlayers = triggerTools_.layers(ForwardSubdetector::HFNose);
0092     else {
0093       nlayers = triggerTools_.layers(ForwardSubdetector::ForwardEmpty);
0094     }
0095     if (layer == 0 || layer > nlayers)
0096       continue;
0097     layers[layer - 1] = true;  //layer 0 doesn't exist, so shift by -1
0098   }
0099   int length = 0;
0100   int maxlength = 0;
0101   for (bool layer : layers) {
0102     if (layer)
0103       length++;
0104     else
0105       length = 0;
0106     if (length > maxlength)
0107       maxlength = length;
0108   }
0109   return maxlength;
0110 }
0111
0112 float HGCalShowerShape::percentileLayer(const l1t::HGCalMulticluster& c3d,
0113                                         const HGCalTriggerGeometryBase& triggerGeometry,
0114                                         float quantile) const {
0115   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0116   unsigned nlayers = triggerTools_.layers(ForwardSubdetector::ForwardEmpty);
0117   std::vector<double> layers(nlayers, 0);
0118   for (const auto& id_clu : clustersPtrs) {
0119     const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& triggerCells = id_clu.second->constituents();
0120
0121     for (const auto& id_tc : triggerCells) {
0122       if (!pass(*id_tc.second, c3d))
0123         continue;
0124       unsigned layer = triggerGeometry.triggerLayer(id_tc.second->detId());
0125       if (triggerTools_.isNose(id_tc.second->detId()))
0126         nlayers = triggerTools_.layers(ForwardSubdetector::HFNose);
0127       else {
0128         nlayers = triggerTools_.layers(ForwardSubdetector::ForwardEmpty);
0129       }
0130       if (layer == 0 || layer > nlayers)
0131         continue;
0132       layers[layer - 1] += id_tc.second->pt();  //layer 0 doesn't exist, so shift by -1
0133     }
0134   }
0135   std::partial_sum(layers.begin(), layers.end(), layers.begin());
0136   double pt_threshold = layers.back() * quantile;
0137   unsigned percentile = 0;
0138   for (double pt : layers) {
0139     if (pt > pt_threshold) {
0140       break;
0141     }
0142     percentile++;
0143   }
0144   // Linear interpolation of percentile value
0145   double pt0 = (percentile > 0 ? layers[percentile - 1] : 0.);
0146   double pt1 = (percentile < layers.size() ? layers[percentile] : layers.back());
0147   return percentile + (pt1 - pt0 > 0. ? (pt_threshold - pt0) / (pt1 - pt0) : 0.);
0148 }
0149
0150 float HGCalShowerShape::percentileTriggerCells(const l1t::HGCalMulticluster& c3d, float quantile) const {
0151   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0152   std::set<double> ordered_tcs;
0153   double pt_sum = 0.;
0154   for (const auto& id_clu : clustersPtrs) {
0155     const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& triggerCells = id_clu.second->constituents();
0156
0157     for (const auto& id_tc : triggerCells) {
0158       if (!pass(*id_tc.second, c3d))
0159         continue;
0160       ordered_tcs.emplace(id_tc.second->pt());
0161       pt_sum += id_tc.second->pt();
0162     }
0163   }
0164   double pt_threshold = pt_sum * quantile;
0165   double partial_sum = 0.;
0166   double partial_sum_prev = 0.;
0167   int ntc = 0;
0168   for (auto itr = ordered_tcs.rbegin(); itr != ordered_tcs.rend(); ++itr) {
0169     partial_sum_prev = partial_sum;
0170     partial_sum += *itr;
0171     ntc++;
0172     if (partial_sum > pt_threshold) {
0173       break;
0174     }
0175   }
0176   // Linear interpolation of ntc
0177   return ntc - 1 +
0178          (partial_sum - partial_sum_prev > 0. ? (pt_threshold - partial_sum_prev) / (partial_sum - partial_sum_prev)
0179                                               : 0.);
0180 }
0181
0182 float HGCalShowerShape::sigmaEtaEtaTot(const l1t::HGCalMulticluster& c3d) const { return sqrt(varEtaEta(c3d)); }
0183
0184 float HGCalShowerShape::varEtaEta(const l1t::HGCalMulticluster& c3d) const {
0185   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0186
0187   std::vector<std::pair<float, float>> tc_energy_eta;
0188
0189   for (const auto& id_clu : clustersPtrs) {
0190     const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& triggerCells = id_clu.second->constituents();
0191
0192     for (const auto& id_tc : triggerCells) {
0193       if (!pass(*id_tc.second, c3d))
0194         continue;
0195       tc_energy_eta.emplace_back(std::make_pair(id_tc.second->energy(), id_tc.second->eta()));
0196     }
0197   }
0198
0199   float varee = varXX(tc_energy_eta, c3d.eta());
0200
0201   return varee;
0202 }
0203
0204 float HGCalShowerShape::sigmaPhiPhiTot(const l1t::HGCalMulticluster& c3d) const { return sqrt(varPhiPhi(c3d)); }
0205
0206 float HGCalShowerShape::varPhiPhi(const l1t::HGCalMulticluster& c3d) const {
0207   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0208
0209   std::vector<std::pair<float, float>> tc_energy_phi;
0210
0211   for (const auto& id_clu : clustersPtrs) {
0212     const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& triggerCells = id_clu.second->constituents();
0213
0214     for (const auto& id_tc : triggerCells) {
0215       if (!pass(*id_tc.second, c3d))
0216         continue;
0217       tc_energy_phi.emplace_back(std::make_pair(id_tc.second->energy(), id_tc.second->phi()));
0218     }
0219   }
0220
0221   float varphiphi = varPhiPhi(tc_energy_phi, c3d.phi());
0222
0223   return varphiphi;
0224 }
0225
0226 float HGCalShowerShape::sigmaRRTot(const l1t::HGCalMulticluster& c3d) const { return sqrt(varRR(c3d)); }
0227
0228 float HGCalShowerShape::varRR(const l1t::HGCalMulticluster& c3d) const {
0229   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0230
0231   std::vector<std::pair<float, float>> tc_energy_r;
0232
0233   for (const auto& id_clu : clustersPtrs) {
0234     const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& triggerCells = id_clu.second->constituents();
0235
0236     for (const auto& id_tc : triggerCells) {
0237       if (!pass(*id_tc.second, c3d))
0238         continue;
0239       float r = (id_tc.second->position().z() != 0.
0240                      ? std::sqrt(pow(id_tc.second->position().x(), 2) + pow(id_tc.second->position().y(), 2)) /
0241                            std::abs(id_tc.second->position().z())
0242                      : 0.);
0243       tc_energy_r.emplace_back(std::make_pair(id_tc.second->energy(), r));
0244     }
0245   }
0246
0247   float r_mean = meanX(tc_energy_r);
0248   float vrr = varXX(tc_energy_r, r_mean);
0249
0250   return vrr;
0251 }
0252
0253 float HGCalShowerShape::sigmaEtaEtaMax(const l1t::HGCalMulticluster& c3d) const {
0254   std::unordered_map<int, std::vector<std::pair<float, float>>> tc_layer_energy_eta;
0255   std::unordered_map<int, LorentzVector> layer_LV;
0256
0257   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0258
0259   for (const auto& id_clu : clustersPtrs) {
0260     unsigned layer = triggerTools_.layerWithOffset(id_clu.second->detId());
0261
0262     layer_LV[layer] += id_clu.second->p4();
0263
0264     const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& triggerCells = id_clu.second->constituents();
0265
0266     for (const auto& id_tc : triggerCells) {
0267       if (!pass(*id_tc.second, c3d))
0268         continue;
0269       tc_layer_energy_eta[layer].emplace_back(std::make_pair(id_tc.second->energy(), id_tc.second->eta()));
0270     }
0271   }
0272
0273   float SigmaEtaEtaMax = 0;
0274
0275   for (auto& tc_iter : tc_layer_energy_eta) {
0276     const std::vector<std::pair<float, float>>& energy_eta_layer = tc_iter.second;
0277     const LorentzVector& LV_layer = layer_LV[tc_iter.first];
0278     float SigmaEtaEtaLayer = sigmaXX(energy_eta_layer, LV_layer.eta());  //RMS wrt layer eta, not wrt c3d eta
0279     if (SigmaEtaEtaLayer > SigmaEtaEtaMax)
0280       SigmaEtaEtaMax = SigmaEtaEtaLayer;
0281   }
0282
0283   return SigmaEtaEtaMax;
0284 }
0285
0286 float HGCalShowerShape::sigmaPhiPhiMax(const l1t::HGCalMulticluster& c3d) const {
0287   std::unordered_map<int, std::vector<std::pair<float, float>>> tc_layer_energy_phi;
0288   std::unordered_map<int, LorentzVector> layer_LV;
0289
0290   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0291
0292   for (const auto& id_clu : clustersPtrs) {
0293     unsigned layer = triggerTools_.layerWithOffset(id_clu.second->detId());
0294
0295     layer_LV[layer] += id_clu.second->p4();
0296
0297     const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& triggerCells = id_clu.second->constituents();
0298
0299     for (const auto& id_tc : triggerCells) {
0300       if (!pass(*id_tc.second, c3d))
0301         continue;
0302       tc_layer_energy_phi[layer].emplace_back(std::make_pair(id_tc.second->energy(), id_tc.second->phi()));
0303     }
0304   }
0305
0306   float SigmaPhiPhiMax = 0;
0307
0308   for (auto& tc_iter : tc_layer_energy_phi) {
0309     const std::vector<std::pair<float, float>>& energy_phi_layer = tc_iter.second;
0310     const LorentzVector& LV_layer = layer_LV[tc_iter.first];
0311     float SigmaPhiPhiLayer = sigmaPhiPhi(energy_phi_layer, LV_layer.phi());  //RMS wrt layer phi, not wrt c3d phi
0312     if (SigmaPhiPhiLayer > SigmaPhiPhiMax)
0313       SigmaPhiPhiMax = SigmaPhiPhiLayer;
0314   }
0315
0316   return SigmaPhiPhiMax;
0317 }
0318
0319 float HGCalShowerShape::sigmaRRMax(const l1t::HGCalMulticluster& c3d) const {
0320   std::unordered_map<int, std::vector<std::pair<float, float>>> tc_layer_energy_r;
0321
0322   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0323
0324   for (const auto& id_clu : clustersPtrs) {
0325     unsigned layer = triggerTools_.layerWithOffset(id_clu.second->detId());
0326
0327     const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& triggerCells = id_clu.second->constituents();
0328
0329     for (const auto& id_tc : triggerCells) {
0330       if (!pass(*id_tc.second, c3d))
0331         continue;
0332       float r = (id_tc.second->position().z() != 0.
0333                      ? std::sqrt(pow(id_tc.second->position().x(), 2) + pow(id_tc.second->position().y(), 2)) /
0334                            std::abs(id_tc.second->position().z())
0335                      : 0.);
0336       tc_layer_energy_r[layer].emplace_back(std::make_pair(id_tc.second->energy(), r));
0337     }
0338   }
0339
0340   float SigmaRRMax = 0;
0341
0342   for (auto& tc_iter : tc_layer_energy_r) {
0343     const std::vector<std::pair<float, float>>& energy_r_layer = tc_iter.second;
0344     float r_mean_layer = meanX(energy_r_layer);
0345     float SigmaRRLayer = sigmaXX(energy_r_layer, r_mean_layer);
0346     if (SigmaRRLayer > SigmaRRMax)
0347       SigmaRRMax = SigmaRRLayer;
0348   }
0349
0350   return SigmaRRMax;
0351 }
0352
0353 float HGCalShowerShape::sigmaRRMean(const l1t::HGCalMulticluster& c3d, float radius) const {
0354   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0355   // group trigger cells by layer
0356   std::unordered_map<int, std::vector<edm::Ptr<l1t::HGCalTriggerCell>>> layers_tcs;
0357   for (const auto& id_clu : clustersPtrs) {
0358     unsigned layer = triggerTools_.layerWithOffset(id_clu.second->detId());
0359     const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& triggerCells = id_clu.second->constituents();
0360     for (const auto& id_tc : triggerCells) {
0361       if (!pass(*id_tc.second, c3d))
0362         continue;
0363       layers_tcs[layer].emplace_back(id_tc.second);
0364     }
0365   }
0366
0367   // Select trigger cells within X cm of the max TC in the layer
0368   std::unordered_map<int, std::vector<std::pair<float, float>>> tc_layers_energy_r;
0369   for (const auto& layer_tcs : layers_tcs) {
0370     int layer = layer_tcs.first;
0371     edm::Ptr<l1t::HGCalTriggerCell> max_tc = layer_tcs.second.front();
0372     for (const auto& tc : layer_tcs.second) {
0373       if (tc->energy() > max_tc->energy())
0374         max_tc = tc;
0375     }
0376     for (const auto& tc : layer_tcs.second) {
0377       double dx = tc->position().x() - max_tc->position().x();
0378       double dy = tc->position().y() - max_tc->position().y();
0379       double distance_to_max = std::sqrt(dx * dx + dy * dy);
0380       if (distance_to_max < radius) {
0381         float r = (tc->position().z() != 0.
0382                        ? std::sqrt(tc->position().x() * tc->position().x() + tc->position().y() * tc->position().y()) /
0383                              std::abs(tc->position().z())
0384                        : 0.);
0385         tc_layers_energy_r[layer].emplace_back(std::make_pair(tc->energy(), r));
0386       }
0387     }
0388   }
0389
0390   // Compute srr layer by layer
0391   std::vector<std::pair<float, float>> layers_energy_srr2;
0392   for (const auto& layer_energy_r : tc_layers_energy_r) {
0393     const auto& energies_r = layer_energy_r.second;
0394     float r_mean_layer = meanX(energies_r);
0395     float srr = sigmaXX(energies_r, r_mean_layer);
0396     double energy_sum = 0.;
0397     for (const auto& energy_r : energies_r) {
0398       energy_sum += energy_r.first;
0399     }
0400     layers_energy_srr2.emplace_back(std::make_pair(energy_sum, srr * srr));
0401   }
0402   // Combine all layer srr
0403   float srr2_mean = meanX(layers_energy_srr2);
0404   return std::sqrt(srr2_mean);
0405 }
0406
0407 float HGCalShowerShape::eMax(const l1t::HGCalMulticluster& c3d) const {
0408   std::unordered_map<int, float> layer_energy;
0409
0410   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0411
0412   for (const auto& id_clu : clustersPtrs) {
0413     if (!pass(*id_clu.second, c3d))
0414       continue;
0415     unsigned layer = triggerTools_.layerWithOffset(id_clu.second->detId());
0416     layer_energy[layer] += id_clu.second->energy();
0417   }
0418
0419   float EMax = 0;
0420
0421   for (const auto& layer : layer_energy) {
0422     if (layer.second > EMax)
0423       EMax = layer.second;
0424   }
0425
0426   return EMax;
0427 }
0428
0429 float HGCalShowerShape::meanZ(const l1t::HGCalMulticluster& c3d) const {
0430   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0431
0432   std::vector<std::pair<float, float>> tc_energy_z;
0433
0434   for (const auto& id_clu : clustersPtrs) {
0435     const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& triggerCells = id_clu.second->constituents();
0436
0437     for (const auto& id_tc : triggerCells) {
0438       if (!pass(*id_tc.second, c3d))
0439         continue;
0440       tc_energy_z.emplace_back(id_tc.second->energy(), id_tc.second->position().z());
0441     }
0442   }
0443
0444   return meanX(tc_energy_z);
0445 }
0446
0447 float HGCalShowerShape::sigmaZZ(const l1t::HGCalMulticluster& c3d) const { return sqrt(varZZ(c3d)); }
0448 float HGCalShowerShape::varZZ(const l1t::HGCalMulticluster& c3d) const {
0449   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0450
0451   std::vector<std::pair<float, float>> tc_energy_z;
0452
0453   for (const auto& id_clu : clustersPtrs) {
0454     const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& triggerCells = id_clu.second->constituents();
0455
0456     for (const auto& id_tc : triggerCells) {
0457       if (!pass(*id_tc.second, c3d))
0458         continue;
0459       tc_energy_z.emplace_back(std::make_pair(id_tc.second->energy(), id_tc.second->position().z()));
0460     }
0461   }
0462
0463   float z_mean = meanX(tc_energy_z);
0464   float varzz = varXX(tc_energy_z, z_mean);
0465
0466   return varzz;
0467 }
0468 float HGCalShowerShape::sigmaEtaEtaTot(const l1t::HGCalCluster& c2d) const {
0469   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& cellsPtrs = c2d.constituents();
0470
0471   std::vector<std::pair<float, float>> tc_energy_eta;
0472
0473   for (const auto& id_cell : cellsPtrs) {
0474     if (!pass(*id_cell.second, c2d))
0475       continue;
0476     tc_energy_eta.emplace_back(std::make_pair(id_cell.second->energy(), id_cell.second->eta()));
0477   }
0478
0479   float See = sigmaXX(tc_energy_eta, c2d.eta());
0480
0481   return See;
0482 }
0483
0484 float HGCalShowerShape::sigmaPhiPhiTot(const l1t::HGCalCluster& c2d) const {
0485   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& cellsPtrs = c2d.constituents();
0486
0487   std::vector<std::pair<float, float>> tc_energy_phi;
0488
0489   for (const auto& id_cell : cellsPtrs) {
0490     if (!pass(*id_cell.second, c2d))
0491       continue;
0492     tc_energy_phi.emplace_back(std::make_pair(id_cell.second->energy(), id_cell.second->phi()));
0493   }
0494
0495   float Spp = sigmaPhiPhi(tc_energy_phi, c2d.phi());
0496
0497   return Spp;
0498 }
0499
0500 float HGCalShowerShape::sigmaRRTot(const l1t::HGCalCluster& c2d) const {
0501   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& cellsPtrs = c2d.constituents();
0502
0503   std::vector<std::pair<float, float>> tc_energy_r;
0504
0505   for (const auto& id_cell : cellsPtrs) {
0506     if (!pass(*id_cell.second, c2d))
0507       continue;
0508     float r = (id_cell.second->position().z() != 0.
0509                    ? std::sqrt(pow(id_cell.second->position().x(), 2) + pow(id_cell.second->position().y(), 2)) /
0510                          std::abs(id_cell.second->position().z())
0511                    : 0.);
0512     tc_energy_r.emplace_back(std::make_pair(id_cell.second->energy(), r));
0513   }
0514
0515   float r_mean = meanX(tc_energy_r);
0516   float Srr = sigmaXX(tc_energy_r, r_mean);
0517
0518   return Srr;
0519 }
0520
0521 float HGCalShowerShape::sumLayers(const l1t::HGCalMulticluster& c3d, int start, int end) const {
0522   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0523   unsigned nlayers = triggerTools_.getTriggerGeometry()->lastTriggerLayer();
0524   std::vector<double> layers(nlayers, 0);
0525   for (const auto& id_clu : clustersPtrs) {
0526     const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& triggerCells = id_clu.second->constituents();
0527
0528     for (const auto& id_tc : triggerCells) {
0529       if (!pass(*id_tc.second, c3d))
0530         continue;
0531       unsigned layer = triggerTools_.triggerLayer(id_tc.second->detId());
0532       if (layer == 0 || layer > nlayers)
0533         continue;
0534       layers[layer - 1] += id_tc.second->pt();  //shift by -1 because layer 0 doesn't exist
0535     }
0536   }
0537   double sum_pt = 0;
0538   for (int i = start - 1; i <= end - 1; i++) {  //shift by -1 because layer 1 is layers[0]
0539     sum_pt += layers[i];
0540   }
0541   double tot = 0;
0542   for (unsigned i = 0; i < layers.size(); ++i) {
0543     tot += layers[i];
0544   }
0545   float frac = 0;
0546   if (tot > 0) {
0547     frac = sum_pt / tot;
0548   }
0549   return frac;
0550 }
0551
0552 int HGCalShowerShape::bitmap(const l1t::HGCalMulticluster& c3d, int start, int end, float threshold) const {
0553   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0554   unsigned nlayers = triggerTools_.getTriggerGeometry()->lastTriggerLayer();
0555   std::vector<double> layers(nlayers, 0);
0556   for (const auto& id_clu : clustersPtrs) {
0557     const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& triggerCells = id_clu.second->constituents();
0558
0559     for (const auto& id_tc : triggerCells) {
0560       if (!pass(*id_tc.second, c3d))
0561         continue;
0562       unsigned layer = triggerTools_.triggerLayer(id_tc.second->detId());
0563       if (layer == 0 || layer > nlayers)
0564         continue;
0565       layers[layer - 1] += id_tc.second->pt();  //shift by -1 because layer 0 doesn't exist
0566     }
0567   }
0568   uint32_t bitmap = 0;
0569   const int bitmap_size = 32;
0570   if (start == 0) {
0571     edm::LogWarning("DataNotFound") << "Trying to read layer 0 that doesn't exist, defaulted start to layer 1";
0572     start = 1;
0573   }
0574   if ((end - start) + 1 > bitmap_size) {
0575     edm::LogWarning("TooMuchData") << " Specified bounds cannot fit into bitmap size, defaulting to 0.";
0576   } else {
0577     for (int i = start; i <= end; i++) {
0578       bitmap += (layers[i - 1] > threshold) << (end - (i));
0579     }
0580   }
0581   return bitmap;
0582 }
0583
0584 void HGCalShowerShape::fillShapes(l1t::HGCalMulticluster& c3d, const HGCalTriggerGeometryBase& triggerGeometry) const {
0585   unsigned hcal_offset = triggerTools_.layers(ForwardSubdetector::HGCEE) / 2;
0586   unsigned lastlayer = triggerGeometry.lastTriggerLayer();
0587   const unsigned showermaxlayer = 7;
0588   c3d.setShowerLength(showerLength(c3d));
0589   c3d.setCoreShowerLength(coreShowerLength(c3d, triggerGeometry));
0590   c3d.setFirstLayer(firstLayer(c3d));
0591   c3d.setMaxLayer(maxLayer(c3d));
0592   c3d.setSigmaEtaEtaTot(sigmaEtaEtaTot(c3d));
0593   c3d.setSigmaEtaEtaMax(sigmaEtaEtaMax(c3d));
0594   c3d.setVarEtaEta(varEtaEta(c3d));
0595   c3d.setSigmaPhiPhiTot(sigmaPhiPhiTot(c3d));
0596   c3d.setSigmaPhiPhiMax(sigmaPhiPhiMax(c3d));
0597   c3d.setVarPhiPhi(varPhiPhi(c3d));
0598   c3d.setSigmaZZ(sigmaZZ(c3d));
0599   c3d.setVarZZ(varZZ(c3d));
0600   c3d.setSigmaRRTot(sigmaRRTot(c3d));
0601   c3d.setVarRR(varRR(c3d));
0602   c3d.setSigmaRRMax(sigmaRRMax(c3d));
0603   c3d.setSigmaRRMean(sigmaRRMean(c3d));
0604   c3d.setEMax(eMax(c3d));
0605   c3d.setZBarycenter(meanZ(c3d));
0606   c3d.setLayer10percent(percentileLayer(c3d, triggerGeometry, 0.10));
0607   c3d.setLayer50percent(percentileLayer(c3d, triggerGeometry, 0.50));
0608   c3d.setLayer90percent(percentileLayer(c3d, triggerGeometry, 0.90));
0609   c3d.setTriggerCells67percent(percentileTriggerCells(c3d, 0.67));
0610   c3d.setTriggerCells90percent(percentileTriggerCells(c3d, 0.90));
0611   c3d.setFirst1layers(sumLayers(c3d, 1, 1));
0612   c3d.setFirst3layers(sumLayers(c3d, 1, 3));
0613   c3d.setFirst5layers(sumLayers(c3d, 1, 5));
0614   c3d.setFirstHcal1layers(sumLayers(c3d, hcal_offset, hcal_offset));
0615   c3d.setFirstHcal3layers(sumLayers(c3d, hcal_offset, hcal_offset + 2));
0616   c3d.setFirstHcal5layers(sumLayers(c3d, hcal_offset, hcal_offset + 4));
0617   c3d.setLast1layers(sumLayers(c3d, lastlayer, lastlayer));
0618   c3d.setLast3layers(sumLayers(c3d, lastlayer - 2, lastlayer));
0619   c3d.setLast5layers(sumLayers(c3d, lastlayer - 4, lastlayer));
0620   c3d.setEmax1layers(sumLayers(c3d, showermaxlayer, showermaxlayer));
0621   c3d.setEmax3layers(sumLayers(c3d, showermaxlayer - 1, showermaxlayer + 1));
0622   c3d.setEmax5layers(sumLayers(c3d, showermaxlayer - 2, showermaxlayer + 2));
0623   c3d.setEot(sumLayers(c3d, 1, hcal_offset));
0624   c3d.setEbm0(bitmap(c3d, 1, hcal_offset, 0));
0625   c3d.setEbm1(bitmap(c3d, 1, hcal_offset, 1));
0626   c3d.setHbm(bitmap(c3d, hcal_offset, lastlayer, 0));
0627 }