File indexing completed on 2024-04-06 12:20:38
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 float varZZ(const l1t::HGCalMulticluster& c3d) const;
0042
0043 float sigmaEtaEtaTot(const l1t::HGCalMulticluster& c3d) const;
0044 float sigmaEtaEtaTot(const l1t::HGCalCluster& c2d) const;
0045 float sigmaEtaEtaMax(const l1t::HGCalMulticluster& c3d) const;
0046 float varEtaEta(const l1t::HGCalMulticluster& c3d) const;
0047
0048 float sigmaPhiPhiTot(const l1t::HGCalMulticluster& c3d) const;
0049 float sigmaPhiPhiTot(const l1t::HGCalCluster& c2d) const;
0050 float sigmaPhiPhiMax(const l1t::HGCalMulticluster& c3d) const;
0051 float varPhiPhi(const l1t::HGCalMulticluster& c3d) const;
0052
0053 float sigmaRRTot(const l1t::HGCalMulticluster& c3d) const;
0054 float sigmaRRTot(const l1t::HGCalCluster& c2d) const;
0055 float sigmaRRMax(const l1t::HGCalMulticluster& c3d) const;
0056 float sigmaRRMean(const l1t::HGCalMulticluster& c3d, float radius = 5.) const;
0057 float varRR(const l1t::HGCalMulticluster& c3d) const;
0058 float sumLayers(const l1t::HGCalMulticluster& c3d, int start = 1, int end = 0) const;
0059 int bitmap(const l1t::HGCalMulticluster& c3d, int start = 1, int end = 14, float threshold = 0) const;
0060 void fillShapes(l1t::HGCalMulticluster&, const HGCalTriggerGeometryBase&) const;
0061
0062 private:
0063 float meanX(const std::vector<pair<float, float>>& energy_X_tc) const;
0064
0065
0066 template <typename Delta = std::minus<float>>
0067 float varXX(const std::vector<pair<float, float>>& energy_X_tc, const float X_cluster) const {
0068 Delta delta;
0069 float Etot = 0;
0070 float deltaX2_sum = 0;
0071 for (const auto& energy_X : energy_X_tc) {
0072 deltaX2_sum += energy_X.first * pow(delta(energy_X.second, X_cluster), 2);
0073 Etot += energy_X.first;
0074 }
0075 float X_MSE = 0;
0076 if (Etot > 0)
0077 X_MSE = deltaX2_sum / Etot;
0078 return X_MSE;
0079 }
0080
0081 float sigmaXX(const std::vector<pair<float, float>>& energy_X_tc, const float X_cluster) const {
0082 float X_RMS = sqrt(varXX(energy_X_tc, X_cluster));
0083 return X_RMS;
0084 }
0085
0086
0087 float sigmaPhiPhi(const std::vector<pair<float, float>>& energy_phi_tc, const float phi_cluster) const {
0088 return sqrt(varPhiPhi(energy_phi_tc, phi_cluster));
0089 }
0090 template <class T>
0091 struct DeltaPhi {
0092 T operator()(const T& x, const T& y) const { return deltaPhi(x, y); }
0093 };
0094 float varPhiPhi(const std::vector<pair<float, float>>& energy_phi_tc, const float phi_cluster) const {
0095 return varXX<DeltaPhi<float>>(energy_phi_tc, phi_cluster);
0096 }
0097
0098 template <typename T, typename Tref>
0099 bool pass(const T& obj, const Tref& ref) const {
0100 bool pass_threshold = (obj.mipPt() > threshold_);
0101 GlobalPoint proj(Basic3DVector<float>(obj.position()) / std::abs(obj.position().z()));
0102 bool pass_distance = ((proj - ref.centreProj()).mag() < distance_);
0103 return pass_threshold && pass_distance;
0104 }
0105
0106 HGCalTriggerTools triggerTools_;
0107 double threshold_ = 0.;
0108 double distance_ = 1.;
0109 };
0110
0111 #endif