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