Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:23:04

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 {
0183   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0184 
0185   std::vector<std::pair<float, float>> tc_energy_eta;
0186 
0187   for (const auto& id_clu : clustersPtrs) {
0188     const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& triggerCells = id_clu.second->constituents();
0189 
0190     for (const auto& id_tc : triggerCells) {
0191       if (!pass(*id_tc.second, c3d))
0192         continue;
0193       tc_energy_eta.emplace_back(std::make_pair(id_tc.second->energy(), id_tc.second->eta()));
0194     }
0195   }
0196 
0197   float SeeTot = sigmaXX(tc_energy_eta, c3d.eta());
0198 
0199   return SeeTot;
0200 }
0201 
0202 float HGCalShowerShape::sigmaPhiPhiTot(const l1t::HGCalMulticluster& c3d) const {
0203   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0204 
0205   std::vector<std::pair<float, float>> tc_energy_phi;
0206 
0207   for (const auto& id_clu : clustersPtrs) {
0208     const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& triggerCells = id_clu.second->constituents();
0209 
0210     for (const auto& id_tc : triggerCells) {
0211       if (!pass(*id_tc.second, c3d))
0212         continue;
0213       tc_energy_phi.emplace_back(std::make_pair(id_tc.second->energy(), id_tc.second->phi()));
0214     }
0215   }
0216 
0217   float SppTot = sigmaPhiPhi(tc_energy_phi, c3d.phi());
0218 
0219   return SppTot;
0220 }
0221 
0222 float HGCalShowerShape::sigmaRRTot(const l1t::HGCalMulticluster& c3d) const {
0223   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0224 
0225   std::vector<std::pair<float, float>> tc_energy_r;
0226 
0227   for (const auto& id_clu : clustersPtrs) {
0228     const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& triggerCells = id_clu.second->constituents();
0229 
0230     for (const auto& id_tc : triggerCells) {
0231       if (!pass(*id_tc.second, c3d))
0232         continue;
0233       float r = (id_tc.second->position().z() != 0.
0234                      ? std::sqrt(pow(id_tc.second->position().x(), 2) + pow(id_tc.second->position().y(), 2)) /
0235                            std::abs(id_tc.second->position().z())
0236                      : 0.);
0237       tc_energy_r.emplace_back(std::make_pair(id_tc.second->energy(), r));
0238     }
0239   }
0240 
0241   float r_mean = meanX(tc_energy_r);
0242   float Szz = sigmaXX(tc_energy_r, r_mean);
0243 
0244   return Szz;
0245 }
0246 
0247 float HGCalShowerShape::sigmaEtaEtaMax(const l1t::HGCalMulticluster& c3d) const {
0248   std::unordered_map<int, std::vector<std::pair<float, float>>> tc_layer_energy_eta;
0249   std::unordered_map<int, LorentzVector> layer_LV;
0250 
0251   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0252 
0253   for (const auto& id_clu : clustersPtrs) {
0254     unsigned layer = triggerTools_.layerWithOffset(id_clu.second->detId());
0255 
0256     layer_LV[layer] += id_clu.second->p4();
0257 
0258     const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& triggerCells = id_clu.second->constituents();
0259 
0260     for (const auto& id_tc : triggerCells) {
0261       if (!pass(*id_tc.second, c3d))
0262         continue;
0263       tc_layer_energy_eta[layer].emplace_back(std::make_pair(id_tc.second->energy(), id_tc.second->eta()));
0264     }
0265   }
0266 
0267   float SigmaEtaEtaMax = 0;
0268 
0269   for (auto& tc_iter : tc_layer_energy_eta) {
0270     const std::vector<std::pair<float, float>>& energy_eta_layer = tc_iter.second;
0271     const LorentzVector& LV_layer = layer_LV[tc_iter.first];
0272     float SigmaEtaEtaLayer = sigmaXX(energy_eta_layer, LV_layer.eta());  //RMS wrt layer eta, not wrt c3d eta
0273     if (SigmaEtaEtaLayer > SigmaEtaEtaMax)
0274       SigmaEtaEtaMax = SigmaEtaEtaLayer;
0275   }
0276 
0277   return SigmaEtaEtaMax;
0278 }
0279 
0280 float HGCalShowerShape::sigmaPhiPhiMax(const l1t::HGCalMulticluster& c3d) const {
0281   std::unordered_map<int, std::vector<std::pair<float, float>>> tc_layer_energy_phi;
0282   std::unordered_map<int, LorentzVector> layer_LV;
0283 
0284   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0285 
0286   for (const auto& id_clu : clustersPtrs) {
0287     unsigned layer = triggerTools_.layerWithOffset(id_clu.second->detId());
0288 
0289     layer_LV[layer] += id_clu.second->p4();
0290 
0291     const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& triggerCells = id_clu.second->constituents();
0292 
0293     for (const auto& id_tc : triggerCells) {
0294       if (!pass(*id_tc.second, c3d))
0295         continue;
0296       tc_layer_energy_phi[layer].emplace_back(std::make_pair(id_tc.second->energy(), id_tc.second->phi()));
0297     }
0298   }
0299 
0300   float SigmaPhiPhiMax = 0;
0301 
0302   for (auto& tc_iter : tc_layer_energy_phi) {
0303     const std::vector<std::pair<float, float>>& energy_phi_layer = tc_iter.second;
0304     const LorentzVector& LV_layer = layer_LV[tc_iter.first];
0305     float SigmaPhiPhiLayer = sigmaPhiPhi(energy_phi_layer, LV_layer.phi());  //RMS wrt layer phi, not wrt c3d phi
0306     if (SigmaPhiPhiLayer > SigmaPhiPhiMax)
0307       SigmaPhiPhiMax = SigmaPhiPhiLayer;
0308   }
0309 
0310   return SigmaPhiPhiMax;
0311 }
0312 
0313 float HGCalShowerShape::sigmaRRMax(const l1t::HGCalMulticluster& c3d) const {
0314   std::unordered_map<int, std::vector<std::pair<float, float>>> tc_layer_energy_r;
0315 
0316   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0317 
0318   for (const auto& id_clu : clustersPtrs) {
0319     unsigned layer = triggerTools_.layerWithOffset(id_clu.second->detId());
0320 
0321     const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& triggerCells = id_clu.second->constituents();
0322 
0323     for (const auto& id_tc : triggerCells) {
0324       if (!pass(*id_tc.second, c3d))
0325         continue;
0326       float r = (id_tc.second->position().z() != 0.
0327                      ? std::sqrt(pow(id_tc.second->position().x(), 2) + pow(id_tc.second->position().y(), 2)) /
0328                            std::abs(id_tc.second->position().z())
0329                      : 0.);
0330       tc_layer_energy_r[layer].emplace_back(std::make_pair(id_tc.second->energy(), r));
0331     }
0332   }
0333 
0334   float SigmaRRMax = 0;
0335 
0336   for (auto& tc_iter : tc_layer_energy_r) {
0337     const std::vector<std::pair<float, float>>& energy_r_layer = tc_iter.second;
0338     float r_mean_layer = meanX(energy_r_layer);
0339     float SigmaRRLayer = sigmaXX(energy_r_layer, r_mean_layer);
0340     if (SigmaRRLayer > SigmaRRMax)
0341       SigmaRRMax = SigmaRRLayer;
0342   }
0343 
0344   return SigmaRRMax;
0345 }
0346 
0347 float HGCalShowerShape::sigmaRRMean(const l1t::HGCalMulticluster& c3d, float radius) const {
0348   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0349   // group trigger cells by layer
0350   std::unordered_map<int, std::vector<edm::Ptr<l1t::HGCalTriggerCell>>> layers_tcs;
0351   for (const auto& id_clu : clustersPtrs) {
0352     unsigned layer = triggerTools_.layerWithOffset(id_clu.second->detId());
0353     const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& triggerCells = id_clu.second->constituents();
0354     for (const auto& id_tc : triggerCells) {
0355       if (!pass(*id_tc.second, c3d))
0356         continue;
0357       layers_tcs[layer].emplace_back(id_tc.second);
0358     }
0359   }
0360 
0361   // Select trigger cells within X cm of the max TC in the layer
0362   std::unordered_map<int, std::vector<std::pair<float, float>>> tc_layers_energy_r;
0363   for (const auto& layer_tcs : layers_tcs) {
0364     int layer = layer_tcs.first;
0365     edm::Ptr<l1t::HGCalTriggerCell> max_tc = layer_tcs.second.front();
0366     for (const auto& tc : layer_tcs.second) {
0367       if (tc->energy() > max_tc->energy())
0368         max_tc = tc;
0369     }
0370     for (const auto& tc : layer_tcs.second) {
0371       double dx = tc->position().x() - max_tc->position().x();
0372       double dy = tc->position().y() - max_tc->position().y();
0373       double distance_to_max = std::sqrt(dx * dx + dy * dy);
0374       if (distance_to_max < radius) {
0375         float r = (tc->position().z() != 0.
0376                        ? std::sqrt(tc->position().x() * tc->position().x() + tc->position().y() * tc->position().y()) /
0377                              std::abs(tc->position().z())
0378                        : 0.);
0379         tc_layers_energy_r[layer].emplace_back(std::make_pair(tc->energy(), r));
0380       }
0381     }
0382   }
0383 
0384   // Compute srr layer by layer
0385   std::vector<std::pair<float, float>> layers_energy_srr2;
0386   for (const auto& layer_energy_r : tc_layers_energy_r) {
0387     const auto& energies_r = layer_energy_r.second;
0388     float r_mean_layer = meanX(energies_r);
0389     float srr = sigmaXX(energies_r, r_mean_layer);
0390     double energy_sum = 0.;
0391     for (const auto& energy_r : energies_r) {
0392       energy_sum += energy_r.first;
0393     }
0394     layers_energy_srr2.emplace_back(std::make_pair(energy_sum, srr * srr));
0395   }
0396   // Combine all layer srr
0397   float srr2_mean = meanX(layers_energy_srr2);
0398   return std::sqrt(srr2_mean);
0399 }
0400 
0401 float HGCalShowerShape::eMax(const l1t::HGCalMulticluster& c3d) const {
0402   std::unordered_map<int, float> layer_energy;
0403 
0404   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0405 
0406   for (const auto& id_clu : clustersPtrs) {
0407     if (!pass(*id_clu.second, c3d))
0408       continue;
0409     unsigned layer = triggerTools_.layerWithOffset(id_clu.second->detId());
0410     layer_energy[layer] += id_clu.second->energy();
0411   }
0412 
0413   float EMax = 0;
0414 
0415   for (const auto& layer : layer_energy) {
0416     if (layer.second > EMax)
0417       EMax = layer.second;
0418   }
0419 
0420   return EMax;
0421 }
0422 
0423 float HGCalShowerShape::meanZ(const l1t::HGCalMulticluster& c3d) const {
0424   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0425 
0426   std::vector<std::pair<float, float>> tc_energy_z;
0427 
0428   for (const auto& id_clu : clustersPtrs) {
0429     const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& triggerCells = id_clu.second->constituents();
0430 
0431     for (const auto& id_tc : triggerCells) {
0432       if (!pass(*id_tc.second, c3d))
0433         continue;
0434       tc_energy_z.emplace_back(id_tc.second->energy(), id_tc.second->position().z());
0435     }
0436   }
0437 
0438   return meanX(tc_energy_z);
0439 }
0440 
0441 float HGCalShowerShape::sigmaZZ(const l1t::HGCalMulticluster& c3d) const {
0442   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalCluster>>& clustersPtrs = c3d.constituents();
0443 
0444   std::vector<std::pair<float, float>> tc_energy_z;
0445 
0446   for (const auto& id_clu : clustersPtrs) {
0447     const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& triggerCells = id_clu.second->constituents();
0448 
0449     for (const auto& id_tc : triggerCells) {
0450       if (!pass(*id_tc.second, c3d))
0451         continue;
0452       tc_energy_z.emplace_back(std::make_pair(id_tc.second->energy(), id_tc.second->position().z()));
0453     }
0454   }
0455 
0456   float z_mean = meanX(tc_energy_z);
0457   float Szz = sigmaXX(tc_energy_z, z_mean);
0458 
0459   return Szz;
0460 }
0461 
0462 float HGCalShowerShape::sigmaEtaEtaTot(const l1t::HGCalCluster& c2d) const {
0463   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& cellsPtrs = c2d.constituents();
0464 
0465   std::vector<std::pair<float, float>> tc_energy_eta;
0466 
0467   for (const auto& id_cell : cellsPtrs) {
0468     if (!pass(*id_cell.second, c2d))
0469       continue;
0470     tc_energy_eta.emplace_back(std::make_pair(id_cell.second->energy(), id_cell.second->eta()));
0471   }
0472 
0473   float See = sigmaXX(tc_energy_eta, c2d.eta());
0474 
0475   return See;
0476 }
0477 
0478 float HGCalShowerShape::sigmaPhiPhiTot(const l1t::HGCalCluster& c2d) const {
0479   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& cellsPtrs = c2d.constituents();
0480 
0481   std::vector<std::pair<float, float>> tc_energy_phi;
0482 
0483   for (const auto& id_cell : cellsPtrs) {
0484     if (!pass(*id_cell.second, c2d))
0485       continue;
0486     tc_energy_phi.emplace_back(std::make_pair(id_cell.second->energy(), id_cell.second->phi()));
0487   }
0488 
0489   float Spp = sigmaPhiPhi(tc_energy_phi, c2d.phi());
0490 
0491   return Spp;
0492 }
0493 
0494 float HGCalShowerShape::sigmaRRTot(const l1t::HGCalCluster& c2d) const {
0495   const std::unordered_map<uint32_t, edm::Ptr<l1t::HGCalTriggerCell>>& cellsPtrs = c2d.constituents();
0496 
0497   std::vector<std::pair<float, float>> tc_energy_r;
0498 
0499   for (const auto& id_cell : cellsPtrs) {
0500     if (!pass(*id_cell.second, c2d))
0501       continue;
0502     float r = (id_cell.second->position().z() != 0.
0503                    ? std::sqrt(pow(id_cell.second->position().x(), 2) + pow(id_cell.second->position().y(), 2)) /
0504                          std::abs(id_cell.second->position().z())
0505                    : 0.);
0506     tc_energy_r.emplace_back(std::make_pair(id_cell.second->energy(), r));
0507   }
0508 
0509   float r_mean = meanX(tc_energy_r);
0510   float Srr = sigmaXX(tc_energy_r, r_mean);
0511 
0512   return Srr;
0513 }
0514 
0515 void HGCalShowerShape::fillShapes(l1t::HGCalMulticluster& c3d, const HGCalTriggerGeometryBase& triggerGeometry) const {
0516   c3d.showerLength(showerLength(c3d));
0517   c3d.coreShowerLength(coreShowerLength(c3d, triggerGeometry));
0518   c3d.firstLayer(firstLayer(c3d));
0519   c3d.maxLayer(maxLayer(c3d));
0520   c3d.sigmaEtaEtaTot(sigmaEtaEtaTot(c3d));
0521   c3d.sigmaEtaEtaMax(sigmaEtaEtaMax(c3d));
0522   c3d.sigmaPhiPhiTot(sigmaPhiPhiTot(c3d));
0523   c3d.sigmaPhiPhiMax(sigmaPhiPhiMax(c3d));
0524   c3d.sigmaZZ(sigmaZZ(c3d));
0525   c3d.sigmaRRTot(sigmaRRTot(c3d));
0526   c3d.sigmaRRMax(sigmaRRMax(c3d));
0527   c3d.sigmaRRMean(sigmaRRMean(c3d));
0528   c3d.eMax(eMax(c3d));
0529   c3d.zBarycenter(meanZ(c3d));
0530   c3d.layer10percent(percentileLayer(c3d, triggerGeometry, 0.10));
0531   c3d.layer50percent(percentileLayer(c3d, triggerGeometry, 0.50));
0532   c3d.layer90percent(percentileLayer(c3d, triggerGeometry, 0.90));
0533   c3d.triggerCells67percent(percentileTriggerCells(c3d, 0.67));
0534   c3d.triggerCells90percent(percentileTriggerCells(c3d, 0.90));
0535 }