interface/backend/HGCalShowerShape.h

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   }  //in number of layers
0029   // Maximum number of consecutive layers in the cluster
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   // Compute energy-weighted RMS of any variable X in the cluster
0065   // Delta(a,b) functor as template argument. Default is (a-b)
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   // Special case of delta for phi
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