Mustache.cc

Go to the documentation of this file.

#include "RecoEcal/EgammaCoreTools/interface/Mustache.h"
#include "TMath.h"
#include "TVector2.h"
#include <cmath>
using namespace std;

namespace reco {  
  namespace MustacheKernel {    
    bool inMustache(const float maxEta, const float maxPhi, 
            const float ClustE, const float ClusEta, 
            const float ClusPhi){
      //bool inMust=false;
      //float eta0 = maxEta;
      //float phi0 = maxPhi;      
      
      constexpr float p00 = -0.107537;
      constexpr float p01 = 0.590969;
      constexpr float p02 = -0.076494;
      constexpr float p10 = -0.0268843;
      constexpr float p11 = 0.147742;
      constexpr float p12 = -0.0191235;
      
      constexpr float w00 = -0.00571429;
      constexpr float w01 = -0.002;
      constexpr float w10 = 0.0135714;
      constexpr float w11 = 0.001;
      
      const float sineta0 = std::sin(maxEta);
      const float eta0xsineta0 = maxEta*sineta0;
      
      
      //2 parabolas (upper and lower) 
      //of the form: y = a*x*x + b
      
      //b comes from a fit to the width
      //and has a slight dependence on E on the upper edge    
      // this only works because of fine tuning :-D
      const float sqrt_log10_clustE = std::sqrt(std::log10(ClustE)+1.1);
      // we need to have this in two steps, so that we don't improperly shift
      // the lower bound!
      float b_upper = w10*eta0xsineta0 + w11 / sqrt_log10_clustE;      
      float b_lower = w00*eta0xsineta0 + w01 / sqrt_log10_clustE; 
      const float midpoint =  0.5*( b_upper + b_lower );
      b_upper -= midpoint;
      b_lower -= midpoint;

      //the curvature comes from a parabolic 
      //fit for many slices in eta given a 
      //slice -0.1 < log10(Et) < 0.1
      const float curv_up=eta0xsineta0*(p00*eta0xsineta0+p01)+p02;
      const float curv_low=eta0xsineta0*(p10*eta0xsineta0+p11)+p12;
      
      //solving for the curviness given the width of this particular point
      const float a_upper=(1/(4*curv_up))-fabs(b_upper);
      const float a_lower = (1/(4*curv_low))-fabs(b_lower);
      
      const double dphi=TVector2::Phi_mpi_pi(ClusPhi-maxPhi);
      const double dphi2 = dphi*dphi;
      // minimum offset is half a crystal width in either direction
      // because science.
      const float upper_cut=( std::max((1./(4.*a_upper)),0.0)*dphi2 +
                  std::max(b_upper,0.0087f) ) + 0.0087;
      const float lower_cut=( std::max((1./(4.*a_lower)),0.0)*dphi2 + 
                  std::min(b_lower,-0.0087f) );
      

      //if(deta < upper_cut && deta > lower_cut) inMust=true;
      
      const float deta=(1-2*(maxEta<0))*(ClusEta-maxEta); // sign flip deta
      return (deta < upper_cut && deta > lower_cut);
    }

    bool inDynamicDPhiWindow(const float seedEta, const float seedPhi,
                 const float ClustE, const float ClusEta,
                 const float ClusPhi) {
      // from Rishi's fits 06 June 2013 in log base 10
      constexpr double yoffsetEB = 7.151e-02;
      constexpr double scaleEB   = 5.656e-01;
      constexpr double xoffsetEB = 2.931e-01;
      constexpr double widthEB   = 2.976e-01;

      constexpr double yoffsetEE_0 =  5.058e-02;
      constexpr double scaleEE_0   =  7.131e-01;
      constexpr double xoffsetEE_0 =  1.668e-02;
      constexpr double widthEE_0   =  4.114e-01;

      constexpr double yoffsetEE_1 = -9.913e-02;
      constexpr double scaleEE_1   =  4.404e+01;
      constexpr double xoffsetEE_1 = -5.326e+00;
      constexpr double widthEE_1   =  1.184e+00;

      constexpr double yoffsetEE_2 = -6.346e-01;
      constexpr double scaleEE_2   =  1.317e+01;
      constexpr double xoffsetEE_2 = -7.037e+00;
      constexpr double widthEE_2   =  2.836e+00;
      
      
      const double absSeedEta = std::abs(seedEta);
      const int etaBin = ( (int)(absSeedEta >= 1.479) + 
               (int)(absSeedEta >= 1.75)  +
               (int)(absSeedEta >= 2.0)     );
      const double logClustEt = std::log10(ClustE/std::cosh(ClusEta));
      const double clusDphi = std::abs(TVector2::Phi_mpi_pi(seedPhi - 
                                ClusPhi));

      double yoffset, scale, xoffset, width, saturation, cutoff, maxdphi;

      switch( etaBin ) {
      case 0: // EB
    yoffset = yoffsetEB;
    scale   = scaleEB;
    xoffset = xoffsetEB;
    width   = 1.0/widthEB;
    saturation = 0.14;
    cutoff     = 0.60;
    break;
      case 1: // 1.479 -> 1.75
    yoffset = yoffsetEE_0;
    scale   = scaleEE_0;
    xoffset = xoffsetEE_0;
    width   = 1.0/widthEE_0;
    saturation = 0.14;
    cutoff     = 0.55;
    break;
      case 2: // 1.75 -> 2.0
    yoffset = yoffsetEE_1;
    scale   = scaleEE_1;
    xoffset = xoffsetEE_1;
    width   = 1.0/widthEE_1;
    saturation = 0.12;
    cutoff     = 0.45;
    break;
      case 3: // 2.0 and up
    yoffset = yoffsetEE_2;
    scale   = scaleEE_2;
    xoffset = xoffsetEE_2;
    width   = 1.0/widthEE_2;
    saturation = 0.12;
    cutoff     = 0.30;
    break;
      default:
    throw cms::Exception("InValidEtaBin")
      << "Calculated invalid eta bin = " << etaBin 
      << " in \"inDynamicDPhiWindow\"" << std::endl;
      }
      
      maxdphi = yoffset+scale/(1+std::exp((logClustEt-xoffset)*width));
      maxdphi = std::min(maxdphi,cutoff);
      maxdphi = std::max(maxdphi,saturation);
      
      return clusDphi < maxdphi;
    }
  }
  
  void Mustache::MustacheID(const reco::SuperCluster& sc, 
                int & nclusters, 
                float & EoutsideMustache) {
    MustacheID(sc.clustersBegin(),sc.clustersEnd(), 
           nclusters, EoutsideMustache);
  }
  
  void Mustache::MustacheID(const CaloClusterPtrVector& clusters, 
                int & nclusters, 
                float & EoutsideMustache) {    
    MustacheID(clusters.begin(),clusters.end(),nclusters,EoutsideMustache);
  }
  
  void Mustache::MustacheID(const std::vector<const CaloCluster*>& clusters, 
                int & nclusters, 
                float & EoutsideMustache) {
    MustacheID(clusters.cbegin(),clusters.cend(),nclusters,EoutsideMustache);
  }

  template<class RandomAccessPtrIterator>
  void Mustache::MustacheID(const RandomAccessPtrIterator& begin, 
                const RandomAccessPtrIterator& end,
                int & nclusters, 
                float & EoutsideMustache) {    
    nclusters = 0;
    EoutsideMustache = 0;
    
    unsigned int ncl = end-begin;
    if(!ncl) return;
    
    //loop over all clusters to find the one with highest energy
    RandomAccessPtrIterator icl = begin;
    RandomAccessPtrIterator clmax = end;
    float emax = 0;
    for( ; icl != end; ++icl){
      const float e = (*icl)->energy();
      if(e > emax){
    emax = e;
    clmax = icl;
      }
    }
    
    if(end == clmax) return;
    
    float eta0 = (*clmax)->eta();
    float phi0 = (*clmax)->phi();
    

    bool inMust = false;
    icl = begin;
    for( ; icl != end; ++icl ){
      inMust=MustacheKernel::inMustache(eta0, phi0, 
                    (*icl)->energy(), 
                    (*icl)->eta(), 
                    (*icl)->phi());
      
      nclusters += (int)!inMust;
      EoutsideMustache += (!inMust)*((*icl)->energy()); 
    }
  }
  
  void Mustache::MustacheClust(const std::vector<CaloCluster>& clusters, 
                   std::vector<unsigned int>& insideMust, 
                   std::vector<unsigned int>& outsideMust){  
    unsigned int ncl = clusters.size();
    if(!ncl) return;
    
    //loop over all clusters to find the one with highest energy
    float emax = 0;
    int imax = -1;
    for(unsigned int i=0; i<ncl; ++i){
      float e = (clusters[i]).energy();
      if(e > emax){
    emax = e;
    imax = i;
      }
    }
    
    if(imax<0) return;
    float eta0 = (clusters[imax]).eta();
    float phi0 = (clusters[imax]).phi();
    
    
    for(unsigned int k=0; k<ncl; k++){
      
      bool inMust=MustacheKernel::inMustache(eta0, phi0, 
                         (clusters[k]).energy(), 
                         (clusters[k]).eta(), 
                         (clusters[k]).phi());
      //return indices of Clusters outside the Mustache
      if (!(inMust)){
    outsideMust.push_back(k);
      }
      else{//return indices of Clusters inside the Mustache
    insideMust.push_back(k);
      }
    }
  }
  
  void Mustache::FillMustacheVar(const std::vector<CaloCluster>& clusters){
    Energy_In_Mustache_=0;
    Energy_Outside_Mustache_=0;
    LowestClusterEInMustache_=0;
    excluded_=0;
    included_=0;
    std::multimap<float, unsigned int>OrderedClust;
    std::vector<unsigned int> insideMust;
    std::vector<unsigned int> outsideMust;
    MustacheClust(clusters, insideMust, outsideMust);
    included_=insideMust.size(); excluded_=outsideMust.size();
    for(unsigned int i=0; i<insideMust.size(); ++i){
      unsigned int index=insideMust[i];
      Energy_In_Mustache_=clusters[index].energy()+Energy_In_Mustache_;
      OrderedClust.insert(make_pair(clusters[index].energy(), index));
    }
    for(unsigned int i=0; i<outsideMust.size(); ++i){
      unsigned int index=outsideMust[i];
      Energy_Outside_Mustache_=clusters[index].energy()+Energy_Outside_Mustache_;
      Et_Outside_Mustache_=clusters[index].energy()*sin(clusters[index].position().theta())
    +Et_Outside_Mustache_;
    }
    std::multimap<float, unsigned int>::iterator it;
    it=OrderedClust.begin();
    unsigned int lowEindex=(*it).second; 
    LowestClusterEInMustache_=clusters[lowEindex].energy();
    
  }
}