EgammaCoreTools/src/Mustache.cc

0001 #include "RecoEcal/EgammaCoreTools/interface/Mustache.h"
0002 #include "TVector2.h"
0003 #include <cmath>
0004 using namespace std;
0005
0006 namespace reco {
0007   namespace MustacheKernel {
0008     bool inMustache(const EcalMustacheSCParameters* params,
0009                     const float maxEta,
0010                     const float maxPhi,
0011                     const float ClustE,
0012                     const float ClusEta,
0013                     const float ClusPhi) {
0014       const auto log10ClustE = std::log10(ClustE);
0015       const auto parabola_params = params->parabolaParameters(log10ClustE, std::abs(ClusEta));
0016       if (!parabola_params) {
0017         return false;
0018       }
0019
0020       const float sineta0 = std::sin(maxEta);
0021       const float eta0xsineta0 = maxEta * sineta0;
0022
0023       //2 parabolas (upper and lower)
0024       //of the form: y = a*x*x + b
0025
0026       //b comes from a fit to the width
0027       //and has a slight dependence on E on the upper edge
0028       // this only works because of fine tuning :-D
0029       const float sqrt_log10_clustE = std::sqrt(log10ClustE + params->sqrtLogClustETuning());
0030       const float b_upper =
0031           parabola_params->w1Up[0] * eta0xsineta0 + parabola_params->w1Up[1] / sqrt_log10_clustE -
0032           0.5 * (parabola_params->w1Up[0] * eta0xsineta0 + parabola_params->w1Up[1] / sqrt_log10_clustE +
0033                  parabola_params->w0Up[0] * eta0xsineta0 + parabola_params->w0Up[1] / sqrt_log10_clustE);
0034       const float b_lower =
0035           parabola_params->w0Low[0] * eta0xsineta0 + parabola_params->w0Low[1] / sqrt_log10_clustE -
0036           0.5 * (parabola_params->w1Low[0] * eta0xsineta0 + parabola_params->w1Low[1] / sqrt_log10_clustE +
0037                  parabola_params->w0Low[0] * eta0xsineta0 + parabola_params->w0Low[1] / sqrt_log10_clustE);
0038
0039       //the curvature comes from a parabolic
0040       //fit for many slices in eta given a
0041       //slice -0.1 < log10(Et) < 0.1
0042       const float curv_up =
0043           eta0xsineta0 * (parabola_params->pUp[0] * eta0xsineta0 + parabola_params->pUp[1]) + parabola_params->pUp[2];
0044       const float curv_low = eta0xsineta0 * (parabola_params->pLow[0] * eta0xsineta0 + parabola_params->pLow[1]) +
0045                              parabola_params->pLow[2];
0046
0047       //solving for the curviness given the width of this particular point
0048       const float a_upper = (1. / (4. * curv_up)) - std::abs(b_upper);
0049       const float a_lower = (1. / (4. * curv_low)) - std::abs(b_lower);
0050
0051       const double dphi = TVector2::Phi_mpi_pi(ClusPhi - maxPhi);
0052       const double dphi2 = dphi * dphi;
0053       // minimum offset is half a crystal width in either direction
0054       // because science.
0055       constexpr float half_crystal_width = 0.0087;
0056       const float upper_cut =
0057           (std::max((1. / (4. * a_upper)), 0.0) * dphi2 + std::max(b_upper, half_crystal_width)) + half_crystal_width;
0058       const float lower_cut = (std::max((1. / (4. * a_lower)), 0.0) * dphi2 + std::min(b_lower, -half_crystal_width));
0059
0060       const float deta = (1 - 2 * (maxEta < 0)) * (ClusEta - maxEta);  // sign flip deta
0061       return (deta < upper_cut && deta > lower_cut);
0062     }
0063
0064     bool inDynamicDPhiWindow(const EcalSCDynamicDPhiParameters* params,
0065                              const float seedEta,
0066                              const float seedPhi,
0067                              const float ClustE,
0068                              const float ClusEta,
0069                              const float ClusPhi) {
0070       const double absSeedEta = std::abs(seedEta);
0071       const double logClustEt = std::log10(ClustE / std::cosh(ClusEta));
0072       const double clusDphi = std::abs(TVector2::Phi_mpi_pi(seedPhi - ClusPhi));
0073
0074       const auto dynamicDPhiParams = params->dynamicDPhiParameters(ClustE, absSeedEta);
0075       if (!dynamicDPhiParams) {
0076         return false;
0077       }
0078
0079       auto maxdphi = dynamicDPhiParams->yoffset +
0080                      dynamicDPhiParams->scale /
0081                          (1. + std::exp((logClustEt - dynamicDPhiParams->xoffset) / dynamicDPhiParams->width));
0082       maxdphi = std::min(maxdphi, dynamicDPhiParams->cutoff);
0083       maxdphi = std::max(maxdphi, dynamicDPhiParams->saturation);
0084
0085       return clusDphi < maxdphi;
0086     }
0087   }  // namespace MustacheKernel
0088
0089   Mustache::Mustache(const EcalMustacheSCParameters* mustache_params) : mustache_params_(mustache_params) {}
0090
0091   void Mustache::MustacheID(const reco::SuperCluster& sc, int& nclusters, float& EoutsideMustache) {
0092     MustacheID(sc.clustersBegin(), sc.clustersEnd(), nclusters, EoutsideMustache);
0093   }
0094
0095   void Mustache::MustacheID(const CaloClusterPtrVector& clusters, int& nclusters, float& EoutsideMustache) {
0096     MustacheID(clusters.begin(), clusters.end(), nclusters, EoutsideMustache);
0097   }
0098
0099   void Mustache::MustacheID(const std::vector<const CaloCluster*>& clusters, int& nclusters, float& EoutsideMustache) {
0100     MustacheID(clusters.cbegin(), clusters.cend(), nclusters, EoutsideMustache);
0101   }
0102
0103   template <class RandomAccessPtrIterator>
0104   void Mustache::MustacheID(const RandomAccessPtrIterator& begin,
0105                             const RandomAccessPtrIterator& end,
0106                             int& nclusters,
0107                             float& EoutsideMustache) {
0108     nclusters = 0;
0109     EoutsideMustache = 0;
0110
0111     unsigned int ncl = end - begin;
0112     if (!ncl)
0113       return;
0114
0115     //loop over all clusters to find the one with highest energy
0116     RandomAccessPtrIterator icl = begin;
0117     RandomAccessPtrIterator clmax = end;
0118     float emax = 0;
0119     for (; icl != end; ++icl) {
0120       const float e = (*icl)->energy();
0121       if (e > emax) {
0122         emax = e;
0123         clmax = icl;
0124       }
0125     }
0126
0127     if (end == clmax)
0128       return;
0129
0130     float eta0 = (*clmax)->eta();
0131     float phi0 = (*clmax)->phi();
0132
0133     bool inMust = false;
0134     icl = begin;
0135     for (; icl != end; ++icl) {
0136       inMust = MustacheKernel::inMustache(mustache_params_, eta0, phi0, (*icl)->energy(), (*icl)->eta(), (*icl)->phi());
0137
0138       nclusters += (int)!inMust;
0139       EoutsideMustache += (!inMust) * ((*icl)->energy());
0140     }
0141   }
0142
0143   void Mustache::MustacheClust(const std::vector<CaloCluster>& clusters,
0144                                std::vector<unsigned int>& insideMust,
0145                                std::vector<unsigned int>& outsideMust) {
0146     unsigned int ncl = clusters.size();
0147     if (!ncl)
0148       return;
0149
0150     //loop over all clusters to find the one with highest energy
0151     float emax = 0;
0152     int imax = -1;
0153     for (unsigned int i = 0; i < ncl; ++i) {
0154       float e = (clusters[i]).energy();
0155       if (e > emax) {
0156         emax = e;
0157         imax = i;
0158       }
0159     }
0160
0161     if (imax < 0)
0162       return;
0163     float eta0 = (clusters[imax]).eta();
0164     float phi0 = (clusters[imax]).phi();
0165
0166     for (unsigned int k = 0; k < ncl; k++) {
0167       bool inMust = MustacheKernel::inMustache(
0168           mustache_params_, eta0, phi0, (clusters[k]).energy(), (clusters[k]).eta(), (clusters[k]).phi());
0169       //return indices of Clusters outside the Mustache
0170       if (!(inMust)) {
0171         outsideMust.push_back(k);
0172       } else {  //return indices of Clusters inside the Mustache
0173         insideMust.push_back(k);
0174       }
0175     }
0176   }
0177
0178   void Mustache::FillMustacheVar(const std::vector<CaloCluster>& clusters) {
0179     Energy_In_Mustache_ = 0;
0180     Energy_Outside_Mustache_ = 0;
0181     LowestClusterEInMustache_ = 0;
0182     excluded_ = 0;
0183     included_ = 0;
0184     std::multimap<float, unsigned int> OrderedClust;
0185     std::vector<unsigned int> insideMust;
0186     std::vector<unsigned int> outsideMust;
0187     MustacheClust(clusters, insideMust, outsideMust);
0188     included_ = insideMust.size();
0189     excluded_ = outsideMust.size();
0190     for (unsigned int i = 0; i < insideMust.size(); ++i) {
0191       unsigned int index = insideMust[i];
0192       Energy_In_Mustache_ = clusters[index].energy() + Energy_In_Mustache_;
0193       OrderedClust.insert(make_pair(clusters[index].energy(), index));
0194     }
0195     for (unsigned int i = 0; i < outsideMust.size(); ++i) {
0196       unsigned int index = outsideMust[i];
0197       Energy_Outside_Mustache_ = clusters[index].energy() + Energy_Outside_Mustache_;
0198       Et_Outside_Mustache_ = clusters[index].energy() * sin(clusters[index].position().theta()) + Et_Outside_Mustache_;
0199     }
0200     std::multimap<float, unsigned int>::iterator it;
0201     it = OrderedClust.begin();
0202     unsigned int lowEindex = (*it).second;
0203     LowestClusterEInMustache_ = clusters[lowEindex].energy();
0204   }
0205 }  // namespace reco