File indexing completed on 2023-03-17 11:12:17
0001 #ifndef __L1Trigger_L1THGCal_HGCalShowerShape_h__
0002 #define __L1Trigger_L1THGCal_HGCalShowerShape_h__
0003 #include <vector>
0004 #include <functional>
0005 #include <cmath>
0006 #include "DataFormats/L1THGCal/interface/HGCalMulticluster.h"
0007 #include "DataFormats/Math/interface/LorentzVector.h"
0008 #include "L1Trigger/L1THGCal/interface/HGCalTriggerGeometryBase.h"
0009 #include "L1Trigger/L1THGCal/interface/HGCalTriggerTools.h"
0010 #include "DataFormats/Math/interface/deltaPhi.h"
0011
0012 class HGCalShowerShape {
0013 public:
0014 typedef math::XYZTLorentzVector LorentzVector;
0015
0016 HGCalShowerShape() : threshold_(0.) {}
0017 HGCalShowerShape(const edm::ParameterSet& conf);
0018
0019 ~HGCalShowerShape() {}
0020
0021 void setGeometry(const HGCalTriggerGeometryBase* const geom) { triggerTools_.setGeometry(geom); }
0022
0023 int firstLayer(const l1t::HGCalMulticluster& c3d) const;
0024 int lastLayer(const l1t::HGCalMulticluster& c3d) const;
0025 int maxLayer(const l1t::HGCalMulticluster& c3d) const;
0026 int showerLength(const l1t::HGCalMulticluster& c3d) const {
0027 return lastLayer(c3d) - firstLayer(c3d) + 1;
0028 }
0029
0030 int coreShowerLength(const l1t::HGCalMulticluster& c3d, const HGCalTriggerGeometryBase& triggerGeometry) const;
0031 float percentileLayer(const l1t::HGCalMulticluster& c3d,
0032 const HGCalTriggerGeometryBase& triggerGeometry,
0033 float quantile = 0.5) const;
0034
0035 float percentileTriggerCells(const l1t::HGCalMulticluster& c3d, float quantile = 0.5) const;
0036
0037 float eMax(const l1t::HGCalMulticluster& c3d) const;
0038
0039 float meanZ(const l1t::HGCalMulticluster& c3d) const;
0040 float sigmaZZ(const l1t::HGCalMulticluster& c3d) const;
0041
0042 float sigmaEtaEtaTot(const l1t::HGCalMulticluster& c3d) const;
0043 float sigmaEtaEtaTot(const l1t::HGCalCluster& c2d) const;
0044 float sigmaEtaEtaMax(const l1t::HGCalMulticluster& c3d) const;
0045
0046 float sigmaPhiPhiTot(const l1t::HGCalMulticluster& c3d) const;
0047 float sigmaPhiPhiTot(const l1t::HGCalCluster& c2d) const;
0048 float sigmaPhiPhiMax(const l1t::HGCalMulticluster& c3d) const;
0049
0050 float sigmaRRTot(const l1t::HGCalMulticluster& c3d) const;
0051 float sigmaRRTot(const l1t::HGCalCluster& c2d) const;
0052 float sigmaRRMax(const l1t::HGCalMulticluster& c3d) const;
0053 float sigmaRRMean(const l1t::HGCalMulticluster& c3d, float radius = 5.) const;
0054
0055 void fillShapes(l1t::HGCalMulticluster&, const HGCalTriggerGeometryBase&) const;
0056
0057 private:
0058 float meanX(const std::vector<pair<float, float>>& energy_X_tc) const;
0059
0060
0061 template <typename Delta = std::minus<float>>
0062 float sigmaXX(const std::vector<pair<float, float>>& energy_X_tc, const float X_cluster) const {
0063 Delta delta;
0064 float Etot = 0;
0065 float deltaX2_sum = 0;
0066 for (const auto& energy_X : energy_X_tc) {
0067 deltaX2_sum += energy_X.first * pow(delta(energy_X.second, X_cluster), 2);
0068 Etot += energy_X.first;
0069 }
0070 float X_MSE = 0;
0071 if (Etot > 0)
0072 X_MSE = deltaX2_sum / Etot;
0073 float X_RMS = sqrt(X_MSE);
0074 return X_RMS;
0075 }
0076
0077 template <class T>
0078 struct DeltaPhi {
0079 T operator()(const T& x, const T& y) const { return deltaPhi(x, y); }
0080 };
0081 float sigmaPhiPhi(const std::vector<pair<float, float>>& energy_phi_tc, const float phi_cluster) const {
0082 return sigmaXX<DeltaPhi<float>>(energy_phi_tc, phi_cluster);
0083 }
0084 template <typename T, typename Tref>
0085 bool pass(const T& obj, const Tref& ref) const {
0086 bool pass_threshold = (obj.mipPt() > threshold_);
0087 GlobalPoint proj(Basic3DVector<float>(obj.position()) / std::abs(obj.position().z()));
0088 bool pass_distance = ((proj - ref.centreProj()).mag() < distance_);
0089 return pass_threshold && pass_distance;
0090 }
0091
0092 HGCalTriggerTools triggerTools_;
0093 double threshold_ = 0.;
0094 double distance_ = 1.;
0095 };
0096
0097 #endif