CMS 3D CMS Logo

Public Member Functions | Private Types | Private Member Functions | Private Attributes

ClusterShapeAlgo Class Reference

#include <ClusterShapeAlgo.h>

List of all members.

Public Member Functions

reco::ClusterShape Calculate (const reco::BasicCluster &passedCluster, const EcalRecHitCollection *hits, const CaloSubdetectorGeometry *geometry, const CaloSubdetectorTopology *topology)
 ClusterShapeAlgo (const edm::ParameterSet &par)
 ClusterShapeAlgo ()

Private Types

enum  { Eta, Phi }

Private Member Functions

double absZernikeMoment (const reco::BasicCluster &passedCluster, int n, int m, double R0=6.6)
double calc_AbsZernikeMoment (const reco::BasicCluster &passedCluster, int n, int m, double R0)
void Calculate_2ndEnergy (const reco::BasicCluster &passedCluster, const EcalRecHitCollection *hits)
void Calculate_BarrelBasketEnergyFraction (const reco::BasicCluster &passedCluster, const EcalRecHitCollection *hits, const int EtaPhi, const CaloSubdetectorGeometry *geometry)
void Calculate_ComplexZernikeMoments (const reco::BasicCluster &passedCluster)
void Calculate_Covariances (const reco::BasicCluster &passedCluster, const EcalRecHitCollection *hits, const CaloSubdetectorGeometry *geometry)
void Calculate_e2x2 ()
void Calculate_e2x5Bottom ()
void Calculate_e2x5Left ()
void Calculate_e2x5Right ()
void Calculate_e2x5Top ()
void Calculate_e3x2 ()
void Calculate_e3x3 ()
void Calculate_e4x4 ()
void Calculate_e5x5 ()
void Calculate_EnergyDepTopology (const reco::BasicCluster &passedCluster, const EcalRecHitCollection *hits, const CaloSubdetectorGeometry *geometry, bool logW=true)
void Calculate_lat (const reco::BasicCluster &passedCluster)
void Calculate_Polynomials (double rho)
void Calculate_TopEnergy (const reco::BasicCluster &passedCluster, const EcalRecHitCollection *hits)
void Create_Map (const EcalRecHitCollection *hits, const CaloSubdetectorTopology *topology)
double f00 (double r)
double f11 (double r)
double f20 (double r)
double f22 (double r)
double f31 (double r)
double f33 (double r)
double f40 (double r)
double f42 (double r)
double f44 (double r)
double f51 (double r)
double f53 (double r)
double f55 (double r)
double factorial (int n) const
double fast_AbsZernikeMoment (const reco::BasicCluster &passedCluster, int n, int m, double R0)

Private Attributes

double A20_
double A42_
double covEtaEta_
double covEtaPhi_
double covPhiPhi_
double e2nd_
DetId e2ndId_
double e2x2_
int e2x2_Diagonal_X_
int e2x2_Diagonal_Y_
double e2x5Bottom_
double e2x5Left_
double e2x5Right_
double e2x5Top_
double e3x2_
double e3x2Ratio_
double e3x3_
double e4x4_
double e5x5_
double eMax_
DetId eMaxId_
std::vector< double > energyBasketFractionEta_
std::vector< double > energyBasketFractionPhi_
std::vector
< EcalClusterEnergyDeposition
energyDistribution_
std::pair< DetId, double > energyMap_ [5][5]
double etaLat_
std::vector< double > fcn_
double lat_
edm::ParameterSet parameterSet_
double phiLat_

Detailed Description

calculates and creates a ClusterShape object

Author:
Michael A. Balazs, UVa
Version:
Id:
ClusterShapeAlgo.h,v 1.17 2010/11/16 15:09:27 argiro Exp

calculates and creates a ClusterShape object

Author:
Michael A. Balazs, UVa
Version:
Id:
SuperClusterShapeAlgo.h,v 1.1 2008/02/11 13:08:05 kkaadze Exp

Definition at line 37 of file ClusterShapeAlgo.h.


Member Enumeration Documentation

anonymous enum [private]
Enumerator:
Eta 
Phi 

Definition at line 110 of file ClusterShapeAlgo.h.

{ Eta, Phi };

Constructor & Destructor Documentation

ClusterShapeAlgo::ClusterShapeAlgo ( const edm::ParameterSet par)

Definition at line 19 of file ClusterShapeAlgo.cc.

                                                             : 
  parameterSet_(par) {}
ClusterShapeAlgo::ClusterShapeAlgo ( ) [inline]

Definition at line 42 of file ClusterShapeAlgo.h.

{ };

Member Function Documentation

double ClusterShapeAlgo::absZernikeMoment ( const reco::BasicCluster passedCluster,
int  n,
int  m,
double  R0 = 6.6 
) [private]

Definition at line 500 of file ClusterShapeAlgo.cc.

References calc_AbsZernikeMoment(), and fast_AbsZernikeMoment().

Referenced by Calculate_ComplexZernikeMoments().

                                                                   {
  // 1. Check if n,m are correctly
  if ((m>n) || ((n-m)%2 != 0) || (n<0) || (m<0)) return -1;

  // 2. Check if n,R0 are within validity Range :
  // n>20 or R0<2.19cm  just makes no sense !
  if ((n>20) || (R0<=2.19)) return -1;
  if (n<=5) return fast_AbsZernikeMoment(passedCluster,n,m,R0);
  else return calc_AbsZernikeMoment(passedCluster,n,m,R0);
}
double ClusterShapeAlgo::calc_AbsZernikeMoment ( const reco::BasicCluster passedCluster,
int  n,
int  m,
double  R0 
) [private]

Definition at line 561 of file ClusterShapeAlgo.cc.

References funct::cos(), energyDistribution_, factorial(), i, m, phi, funct::pow(), csvReporter::r, query::result, asciidump::s, funct::sin(), and mathSSE::sqrt().

Referenced by absZernikeMoment().

                                                                        {
  double r,ph,e,Re=0,Im=0,f_nm,result;
  double TotalEnergy = passedCluster.energy();
  std::vector< std::pair<DetId, float> > clusterDetIds = passedCluster.hitsAndFractions();
  int clusterSize=energyDistribution_.size();
  if(clusterSize<3) return 0.0;

  for (int i=0; i<clusterSize; i++)
    { 
      r = energyDistribution_[i].r / R0;
      if (r<1) {
        ph = (energyDistribution_[i]).phi;
        e = energyDistribution_[i].deposited_energy;
        f_nm=0;
        for (int s=0; s<=(n-m)/2; s++) {
          if (s%2==0)
            { 
              f_nm = f_nm + factorial(n-s)*pow(r,(double) (n-2*s))/(factorial(s)*factorial((n+m)/2-s)*factorial((n-m)/2-s));
            }else {
              f_nm = f_nm - factorial(n-s)*pow(r,(double) (n-2*s))/(factorial(s)*factorial((n+m)/2-s)*factorial((n-m)/2-s));
            }
        }
        Re = Re + e/TotalEnergy * f_nm * cos( (double) m*ph);
        Im = Im - e/TotalEnergy * f_nm * sin( (double) m*ph);
      }
    }
  result = sqrt(Re*Re+Im*Im);

  return result;
}
reco::ClusterShape ClusterShapeAlgo::Calculate ( const reco::BasicCluster passedCluster,
const EcalRecHitCollection hits,
const CaloSubdetectorGeometry geometry,
const CaloSubdetectorTopology topology 
)

Definition at line 22 of file ClusterShapeAlgo.cc.

References A20_, A42_, Calculate_2ndEnergy(), Calculate_BarrelBasketEnergyFraction(), Calculate_ComplexZernikeMoments(), Calculate_Covariances(), Calculate_e2x2(), Calculate_e2x5Bottom(), Calculate_e2x5Left(), Calculate_e2x5Right(), Calculate_e2x5Top(), Calculate_e3x2(), Calculate_e3x3(), Calculate_e4x4(), Calculate_e5x5(), Calculate_EnergyDepTopology(), Calculate_lat(), Calculate_TopEnergy(), covEtaEta_, covEtaPhi_, covPhiPhi_, Create_Map(), e2nd_, e2ndId_, e2x2_, e2x5Bottom_, e2x5Left_, e2x5Right_, e2x5Top_, e3x2_, e3x2Ratio_, e3x3_, e4x4_, e5x5_, eMax_, eMaxId_, energyBasketFractionEta_, energyBasketFractionPhi_, Eta, etaLat_, lat_, Phi, and phiLat_.

Referenced by CosmicClusterProducer::clusterizeECALPart(), IslandClusterProducer::clusterizeECALPart(), and Pi0FixedMassWindowCalibration::duringLoop().

void ClusterShapeAlgo::Calculate_2ndEnergy ( const reco::BasicCluster passedCluster,
const EcalRecHitCollection hits 
) [private]

Definition at line 82 of file ClusterShapeAlgo.cc.

References e2nd_, e2ndId_, eMaxId_, edm::SortedCollection< T, SORT >::end(), CaloRecHit::energy(), edm::SortedCollection< T, SORT >::find(), and EcalRecHit::id().

Referenced by Calculate().

{
  double e2nd=0;
  DetId e2ndId(0);

  std::vector< std::pair<DetId, float> > clusterDetIds = passedCluster.hitsAndFractions();
  std::vector< std::pair<DetId, float> >::iterator posCurrent;

  EcalRecHit testEcalRecHit;

  for(posCurrent = clusterDetIds.begin(); posCurrent != clusterDetIds.end(); posCurrent++)
  { 
    if (( (*posCurrent).first != DetId(0)) && (hits->find( (*posCurrent).first ) != hits->end()))
    {
      EcalRecHitCollection::const_iterator itt = hits->find( (*posCurrent).first );
      testEcalRecHit = *itt;

      if(testEcalRecHit.energy() * (*posCurrent).second > e2nd && testEcalRecHit.id() != eMaxId_)
      {
        e2nd = testEcalRecHit.energy() * (*posCurrent).second;
        e2ndId = testEcalRecHit.id();
      }
    }
  }

  e2nd_ = e2nd;
  e2ndId_ = e2ndId;
}
void ClusterShapeAlgo::Calculate_BarrelBasketEnergyFraction ( const reco::BasicCluster passedCluster,
const EcalRecHitCollection hits,
const int  EtaPhi,
const CaloSubdetectorGeometry geometry 
) [private]

Definition at line 384 of file ClusterShapeAlgo.cc.

References abs, EcalBarrel, relval_parameters_module::energy, energyBasketFractionEta_, energyBasketFractionPhi_, Eta, edm::SortedCollection< T, SORT >::find(), first, EcalBarrelGeometry::getBasketSizeInPhi(), EcalBarrelGeometry::getEtaBaskets(), i, EBDetId::ieta(), EBDetId::iphi(), EBDetId::MAX_IPHI, EBDetId::MIN_IPHI, and Phi.

Referenced by Calculate().

{
  if(  (hits!=0) && ( ((*hits)[0]).id().subdetId() != EcalBarrel )  ) {
     //std::cout << "No basket correction for endacap!" << std::endl;
     return;
  }

  std::map<int,double> indexedBasketEnergy;
  std::vector< std::pair<DetId, float> > clusterDetIds = passedCluster.hitsAndFractions();
  const EcalBarrelGeometry* subDetGeometry = (const EcalBarrelGeometry*) geometry;

  for(std::vector< std::pair<DetId, float> >::iterator posCurrent = clusterDetIds.begin(); posCurrent != clusterDetIds.end(); posCurrent++)
  {
    int basketIndex = 999;

    if(EtaPhi == Eta) {
      int unsignedIEta = abs(EBDetId( (*posCurrent).first ).ieta());
      std::vector<int> etaBasketSize = subDetGeometry->getEtaBaskets();

      for(unsigned int i = 0; i < etaBasketSize.size(); i++) {
        unsignedIEta -= etaBasketSize[i];
        if(unsignedIEta - 1 < 0)
        {
          basketIndex = i;
          break;
        }
      }
      basketIndex = (basketIndex+1)*(EBDetId( (*posCurrent).first ).ieta() > 0 ? 1 : -1);

    } else if(EtaPhi == Phi) {
      int halfNumBasketInPhi = (EBDetId::MAX_IPHI - EBDetId::MIN_IPHI + 1)/subDetGeometry->getBasketSizeInPhi()/2;

      basketIndex = (EBDetId( (*posCurrent).first ).iphi() - 1)/subDetGeometry->getBasketSizeInPhi()
                  - (EBDetId( (clusterDetIds[0]).first ).iphi() - 1)/subDetGeometry->getBasketSizeInPhi();

      if(basketIndex >= halfNumBasketInPhi)             basketIndex -= 2*halfNumBasketInPhi;
      else if(basketIndex <  -1*halfNumBasketInPhi)     basketIndex += 2*halfNumBasketInPhi;

    } else throw(std::runtime_error("\n\nOh No! Calculate_BarrelBasketEnergyFraction called on invalid index.\n\n"));

    indexedBasketEnergy[basketIndex] += (hits->find( (*posCurrent).first ))->energy();
  }

  std::vector<double> energyFraction;
  for(std::map<int,double>::iterator posCurrent = indexedBasketEnergy.begin(); posCurrent != indexedBasketEnergy.end(); posCurrent++)
  {
    energyFraction.push_back(indexedBasketEnergy[posCurrent->first]/passedCluster.energy());
  }

  switch(EtaPhi)
  {
    case Eta: energyBasketFractionEta_ = energyFraction; break;
    case Phi: energyBasketFractionPhi_ = energyFraction; break;
  }

}
void ClusterShapeAlgo::Calculate_ComplexZernikeMoments ( const reco::BasicCluster passedCluster) [private]

Definition at line 493 of file ClusterShapeAlgo.cc.

References A20_, A42_, and absZernikeMoment().

Referenced by Calculate().

                                                                                            {
  // Calculate only the moments which go into the default cluster shape
  // (moments with m>=2 are the only sensitive to azimuthal shape)
  A20_ = absZernikeMoment(passedCluster,2,0);
  A42_ = absZernikeMoment(passedCluster,4,2);
}
void ClusterShapeAlgo::Calculate_Covariances ( const reco::BasicCluster passedCluster,
const EcalRecHitCollection hits,
const CaloSubdetectorGeometry geometry 
) [private]

Definition at line 312 of file ClusterShapeAlgo.cc.

References covEtaEta_, covEtaPhi_, covPhiPhi_, RecoTauValidation_cfi::denominator, dPhi(), e5x5_, energyMap_, PV3DBase< T, PVType, FrameType >::eta(), CaloSubdetectorGeometry::getGeometry(), edm::ParameterSet::getParameter(), CaloCellGeometry::getPosition(), i, j, funct::log(), max(), parameterSet_, PV3DBase< T, PVType, FrameType >::phi(), Geom::pi(), position, edm::second(), Geom::twoPi(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by Calculate().

{
  if (e5x5_ > 0.)
    {
      double w0_ = parameterSet_.getParameter<double>("W0");
      
      
      // first find energy-weighted mean position - doing it when filling the energy map might save time
      math::XYZVector meanPosition(0.0, 0.0, 0.0);
      for (int i = 0; i < 5; ++i)
        {
          for (int j = 0; j < 5; ++j)
            {
              DetId id = energyMap_[i][j].first;
              if (id != DetId(0))
                {
                  GlobalPoint positionGP = geometry->getGeometry(id)->getPosition();
                  math::XYZVector position(positionGP.x(),positionGP.y(),positionGP.z());
                  meanPosition = meanPosition + energyMap_[i][j].second * position;
                }
            }
        }
      
      meanPosition /= e5x5_;
      
      // now we can calculate the covariances
      double numeratorEtaEta = 0;
      double numeratorEtaPhi = 0;
      double numeratorPhiPhi = 0;
      double denominator     = 0;
      
      for (int i = 0; i < 5; ++i)
        {
          for (int j = 0; j < 5; ++j)
            {
              DetId id = energyMap_[i][j].first;
              if (id != DetId(0))
                {
                  GlobalPoint position = geometry->getGeometry(id)->getPosition();
                  
                  double dPhi = position.phi() - meanPosition.phi();
                  if (dPhi > + Geom::pi()) { dPhi = Geom::twoPi() - dPhi; }
                  if (dPhi < - Geom::pi()) { dPhi = Geom::twoPi() + dPhi; }
                  
                  double dEta = position.eta() - meanPosition.eta();
                  double w = 0.;
                  if ( energyMap_[i][j].second > 0.)
                    w = std::max(0.0, w0_ + log( energyMap_[i][j].second / e5x5_));
                  
                  denominator += w;
                  numeratorEtaEta += w * dEta * dEta;
                  numeratorEtaPhi += w * dEta * dPhi;
                  numeratorPhiPhi += w * dPhi * dPhi;
                }
            }
        }
      
      covEtaEta_ = numeratorEtaEta / denominator;
      covEtaPhi_ = numeratorEtaPhi / denominator;
      covPhiPhi_ = numeratorPhiPhi / denominator;
    }
  else 
    {
      // Warn the user if there was no energy in the cells and return zeroes.
      //       std::cout << "\ClusterShapeAlgo::Calculate_Covariances:  no energy in supplied cells.\n";
      covEtaEta_ = 0;
      covEtaPhi_ = 0;
      covPhiPhi_ = 0;
    }
}
void ClusterShapeAlgo::Calculate_e2x2 ( ) [private]

Definition at line 133 of file ClusterShapeAlgo.cc.

References e2x2_, e2x2_Diagonal_X_, e2x2_Diagonal_Y_, and energyMap_.

Referenced by Calculate().

{
  double e2x2Max = 0;
  double e2x2Test = 0;

  int deltaX=0, deltaY=0;

  for(int corner = 0; corner < 4; corner++)
  {
    switch(corner)
    {
      case 0: deltaX = -1; deltaY = -1; break;
      case 1: deltaX = -1; deltaY =  1; break;
      case 2: deltaX =  1; deltaY = -1; break;
      case 3: deltaX =  1; deltaY =  1; break;
    }

    e2x2Test  = energyMap_[2][2].second;
    e2x2Test += energyMap_[2+deltaY][2].second;
    e2x2Test += energyMap_[2][2+deltaX].second;
    e2x2Test += energyMap_[2+deltaY][2+deltaX].second;

    if(e2x2Test > e2x2Max)
    {
                        e2x2Max = e2x2Test;
                        e2x2_Diagonal_X_ = 2+deltaX;
                        e2x2_Diagonal_Y_ = 2+deltaY;
    }
  }

  e2x2_ = e2x2Max;

}
void ClusterShapeAlgo::Calculate_e2x5Bottom ( ) [private]

Definition at line 288 of file ClusterShapeAlgo.cc.

References e2x5Bottom_, energyMap_, i, and j.

Referenced by Calculate().

{
double e2x5B=0.0;
  for(int i = 0; i <= 4; i++){
    for(int j = 0; j <= 4; j++){

      if(i>2){e2x5B +=energyMap_[i][j].second;}
    }
  }
  e2x5Bottom_=e2x5B;
}
void ClusterShapeAlgo::Calculate_e2x5Left ( ) [private]

Definition at line 277 of file ClusterShapeAlgo.cc.

References e2x5Left_, energyMap_, i, and j.

Referenced by Calculate().

{
double e2x5L=0.0;
  for(int i = 0; i <= 4; i++){
    for(int j = 0; j <= 4; j++){
      if(j<2){e2x5L +=energyMap_[i][j].second;}
    }
  }
  e2x5Left_=e2x5L;
}
void ClusterShapeAlgo::Calculate_e2x5Right ( ) [private]

Definition at line 266 of file ClusterShapeAlgo.cc.

References e2x5Right_, energyMap_, i, and j.

Referenced by Calculate().

{
double e2x5R=0.0;
  for(int i = 0; i <= 4; i++){
    for(int j = 0; j <= 4; j++){
      if(j>2){e2x5R +=energyMap_[i][j].second;}
    }
  }
  e2x5Right_=e2x5R;
}
void ClusterShapeAlgo::Calculate_e2x5Top ( ) [private]

Definition at line 300 of file ClusterShapeAlgo.cc.

References e2x5Top_, energyMap_, i, and j.

Referenced by Calculate().

{
double e2x5T=0.0;
  for(int i = 0; i <= 4; i++){
    for(int j = 0; j <= 4; j++){

      if(i<2){e2x5T +=energyMap_[i][j].second;}
    }
  }
  e2x5Top_=e2x5T;
}
void ClusterShapeAlgo::Calculate_e3x2 ( ) [private]

Definition at line 167 of file ClusterShapeAlgo.cc.

References e3x2_, e3x2Ratio_, energyMap_, and edm::second().

Referenced by Calculate().

{
  double e3x2 = 0.0;
  double e3x2Ratio = 0.0, e3x2RatioNumerator = 0.0, e3x2RatioDenominator = 0.0;

  int e2ndX = 2, e2ndY = 2;
  int deltaY = 0, deltaX = 0;

  double nextEnergy = -999;
  int nextEneryDirection = -1;

  for(int cardinalDirection = 0; cardinalDirection < 4; cardinalDirection++)
  {
    switch(cardinalDirection)
    {
      case 0: deltaX = -1; deltaY =  0; break;
      case 1: deltaX =  1; deltaY =  0; break;
      case 2: deltaX =  0; deltaY = -1; break;
      case 3: deltaX =  0; deltaY =  1; break;
    }
   
    if(energyMap_[2+deltaY][2+deltaX].second >= nextEnergy)
    {
        nextEnergy = energyMap_[2+deltaY][2+deltaX].second;
        nextEneryDirection = cardinalDirection;
       
        e2ndX = 2+deltaX;
        e2ndY = 2+deltaY;
    }
  }
 
  switch(nextEneryDirection)
  {
    case 0: ;
    case 1: deltaX = 0; deltaY = 1; break;
    case 2: ;
    case 3: deltaX = 1; deltaY = 0; break;
  }

  for(int sign = -1; sign <= 1; sign++)
      e3x2 += (energyMap_[2+deltaY*sign][2+deltaX*sign].second + energyMap_[e2ndY+deltaY*sign][e2ndX+deltaX*sign].second);
 
  e3x2RatioNumerator   = energyMap_[e2ndY+deltaY][e2ndX+deltaX].second + energyMap_[e2ndY-deltaY][e2ndX-deltaX].second;
  e3x2RatioDenominator = 0.5 + energyMap_[2+deltaY][2+deltaX].second + energyMap_[2-deltaY][2-deltaX].second;
  e3x2Ratio = e3x2RatioNumerator / e3x2RatioDenominator;

  e3x2_ = e3x2;
  e3x2Ratio_ = e3x2Ratio;
}  
void ClusterShapeAlgo::Calculate_e3x3 ( ) [private]

Definition at line 217 of file ClusterShapeAlgo.cc.

References e3x3_, energyMap_, i, j, and edm::second().

Referenced by Calculate().

{
  double e3x3=0;

  for(int i = 1; i <= 3; i++)
    for(int j = 1; j <= 3; j++)
      e3x3 += energyMap_[j][i].second;

  e3x3_ = e3x3;

}
void ClusterShapeAlgo::Calculate_e4x4 ( ) [private]

Definition at line 229 of file ClusterShapeAlgo.cc.

References e2x2_Diagonal_X_, e2x2_Diagonal_Y_, e4x4_, energyMap_, i, j, and edm::second().

Referenced by Calculate().

{
  double e4x4=0;
        
  int startX=-1, startY=-1;

        switch(e2x2_Diagonal_X_)
        {
                case 1: startX = 0; break;
                case 3: startX = 1; break;
        }

        switch(e2x2_Diagonal_Y_)
        {
                case 1: startY = 0; break;
                case 3: startY = 1; break;
        }

  for(int i = startX; i <= startX+3; i++)
          for(int j = startY; j <= startY+3; j++)
          e4x4 += energyMap_[j][i].second;

  e4x4_ = e4x4;
}
void ClusterShapeAlgo::Calculate_e5x5 ( ) [private]

Definition at line 254 of file ClusterShapeAlgo.cc.

References e5x5_, energyMap_, i, j, and edm::second().

Referenced by Calculate().

{
  double e5x5=0;

  for(int i = 0; i <= 4; i++)
    for(int j = 0; j <= 4; j++)
      e5x5 += energyMap_[j][i].second;

  e5x5_ = e5x5;

}
void ClusterShapeAlgo::Calculate_EnergyDepTopology ( const reco::BasicCluster passedCluster,
const EcalRecHitCollection hits,
const CaloSubdetectorGeometry geometry,
bool  logW = true 
) [private]

Definition at line 593 of file ClusterShapeAlgo.cc.

References EcalClusterEnergyDeposition::deposited_energy, diffTreeTool::diff, edm::SortedCollection< T, SORT >::end(), CaloRecHit::energy(), energyDistribution_, edm::SortedCollection< T, SORT >::find(), CaloSubdetectorGeometry::getGeometry(), edm::ParameterSet::getParameter(), CaloCellGeometry::getPosition(), funct::log(), LogDebug, M_PI, max(), parameterSet_, EcalClusterEnergyDeposition::phi, EcalClusterEnergyDeposition::r, CommonMethods::weight(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by Calculate().

                                                               {
  // resets the energy distribution
  energyDistribution_.clear();

  // init a map of the energy deposition centered on the
  // cluster centroid. This is for momenta calculation only.
  CLHEP::Hep3Vector clVect(passedCluster.position().x(),
                           passedCluster.position().y(),
                           passedCluster.position().z());
  CLHEP::Hep3Vector clDir(clVect);
  clDir*=1.0/clDir.mag();
  // in the transverse plane, axis perpendicular to clusterDir
  CLHEP::Hep3Vector theta_axis(clDir.y(),-clDir.x(),0.0);
  theta_axis *= 1.0/theta_axis.mag();
  CLHEP::Hep3Vector phi_axis = theta_axis.cross(clDir);

  std::vector< std::pair<DetId, float> > clusterDetIds = passedCluster.hitsAndFractions();

  EcalClusterEnergyDeposition clEdep;
  EcalRecHit testEcalRecHit;
  std::vector< std::pair<DetId, float> >::iterator posCurrent;
  // loop over crystals
  for(posCurrent=clusterDetIds.begin(); posCurrent!=clusterDetIds.end(); ++posCurrent) {
    EcalRecHitCollection::const_iterator itt = hits->find( (*posCurrent).first );
    testEcalRecHit=*itt;

    if(( (*posCurrent).first != DetId(0)) && (hits->find( (*posCurrent).first ) != hits->end())) {
      clEdep.deposited_energy = testEcalRecHit.energy();

      // if logarithmic weight is requested, apply cut on minimum energy of the recHit
      if(logW) {
        double w0_ = parameterSet_.getParameter<double>("W0");

        if ( clEdep.deposited_energy == 0 ) {
          LogDebug("ClusterShapeAlgo") << "Crystal has zero energy; skipping... ";
          continue;
        }
        double weight = std::max(0.0, w0_ + log(fabs(clEdep.deposited_energy)/passedCluster.energy()) );
        if(weight==0) {
          LogDebug("ClusterShapeAlgo") << "Crystal has insufficient energy: E = " 
                                       << clEdep.deposited_energy << " GeV; skipping... ";
          continue;
        }
        else LogDebug("ClusterShapeAlgo") << "===> got crystal. Energy = " << clEdep.deposited_energy << " GeV. ";
      }
      DetId id_ = (*posCurrent).first;
      const CaloCellGeometry *this_cell = geometry->getGeometry(id_);
      GlobalPoint cellPos = this_cell->getPosition();
      CLHEP::Hep3Vector gblPos (cellPos.x(),cellPos.y(),cellPos.z()); //surface position?
      // Evaluate the distance from the cluster centroid
      CLHEP::Hep3Vector diff = gblPos - clVect;
      // Important: for the moment calculation, only the "lateral distance" is important
      // "lateral distance" r_i = distance of the digi position from the axis Origin-Cluster Center
      // ---> subtract the projection on clDir
      CLHEP::Hep3Vector DigiVect = diff - diff.dot(clDir)*clDir;
      clEdep.r = DigiVect.mag();
      LogDebug("ClusterShapeAlgo") << "E = " << clEdep.deposited_energy
                                   << "\tdiff = " << diff.mag()
                                   << "\tr = " << clEdep.r;
      clEdep.phi = DigiVect.angle(theta_axis);
      if(DigiVect.dot(phi_axis)<0) clEdep.phi = 2*M_PI - clEdep.phi;
      energyDistribution_.push_back(clEdep);
    }
  } 
}
void ClusterShapeAlgo::Calculate_lat ( const reco::BasicCluster passedCluster) [private]

Definition at line 444 of file ClusterShapeAlgo.cc.

References funct::cos(), energyDistribution_, etaLat_, i, lat_, n, phi, phiLat_, csvReporter::r, funct::sin(), and tmp.

Referenced by Calculate().

                                                                          {

  double r,redmoment=0;
  double phiRedmoment = 0 ;
  double etaRedmoment = 0 ;
  int n,n1,n2,tmp;
  int clusterSize=energyDistribution_.size();
  if (clusterSize<3) {
    etaLat_ = 0.0 ; 
    lat_ = 0.0;
    return; 
  }
  
  n1=0; n2=1;
  if (energyDistribution_[1].deposited_energy > 
      energyDistribution_[0].deposited_energy) 
    {
      tmp=n2; n2=n1; n1=tmp;
    }
  for (int i=2; i<clusterSize; i++) {
    n=i;
    if (energyDistribution_[i].deposited_energy > 
        energyDistribution_[n1].deposited_energy) 
      {
        tmp = n2;
        n2 = n1; n1 = i; n=tmp;
      } else {
        if (energyDistribution_[i].deposited_energy > 
            energyDistribution_[n2].deposited_energy) 
          {
            tmp=n2; n2=i; n=tmp;
          }
      }

    r = energyDistribution_[n].r;
    redmoment += r*r* energyDistribution_[n].deposited_energy;
    double rphi = r * cos (energyDistribution_[n].phi) ;
    phiRedmoment += rphi * rphi * energyDistribution_[n].deposited_energy;
    double reta = r * sin (energyDistribution_[n].phi) ;
    etaRedmoment += reta * reta * energyDistribution_[n].deposited_energy;
  } 
  double e1 = energyDistribution_[n1].deposited_energy;
  double e2 = energyDistribution_[n2].deposited_energy;
  
  lat_ = redmoment/(redmoment+2.19*2.19*(e1+e2));
  phiLat_ = phiRedmoment/(phiRedmoment+2.19*2.19*(e1+e2));
  etaLat_ = etaRedmoment/(etaRedmoment+2.19*2.19*(e1+e2));
}
void ClusterShapeAlgo::Calculate_Polynomials ( double  rho) [private]

Definition at line 662 of file ClusterShapeAlgo.cc.

References f00(), f11(), f20(), f22(), f31(), f33(), f40(), f42(), f44(), f51(), f53(), f55(), and fcn_.

Referenced by fast_AbsZernikeMoment().

                                                       {
  fcn_.push_back(f00(rho));
  fcn_.push_back(f11(rho));
  fcn_.push_back(f20(rho));
  fcn_.push_back(f31(rho));
  fcn_.push_back(f22(rho));
  fcn_.push_back(f33(rho));
  fcn_.push_back(f40(rho));
  fcn_.push_back(f51(rho));
  fcn_.push_back(f42(rho));
  fcn_.push_back(f53(rho));
  fcn_.push_back(f44(rho));
  fcn_.push_back(f55(rho));
}
void ClusterShapeAlgo::Calculate_TopEnergy ( const reco::BasicCluster passedCluster,
const EcalRecHitCollection hits 
) [private]

Definition at line 53 of file ClusterShapeAlgo.cc.

References jptDQMConfig_cff::eMax, eMax_, eMaxId_, edm::SortedCollection< T, SORT >::end(), CaloRecHit::energy(), edm::SortedCollection< T, SORT >::find(), and EcalRecHit::id().

Referenced by Calculate().

{
  double eMax=0;
  DetId eMaxId(0);

  std::vector< std::pair<DetId, float> > clusterDetIds = passedCluster.hitsAndFractions();
  std::vector< std::pair<DetId, float> >::iterator posCurrent;

  EcalRecHit testEcalRecHit;

  for(posCurrent = clusterDetIds.begin(); posCurrent != clusterDetIds.end(); posCurrent++)
  {
    if (((*posCurrent).first != DetId(0)) && (hits->find((*posCurrent).first) != hits->end()))
    {
      EcalRecHitCollection::const_iterator itt = hits->find((*posCurrent).first);
      testEcalRecHit = *itt;

      if(testEcalRecHit.energy() * (*posCurrent).second > eMax)
      {
        eMax = testEcalRecHit.energy() * (*posCurrent).second;
        eMaxId = testEcalRecHit.id();
      }
    }
  }

  eMax_ = eMax;
  eMaxId_ = eMaxId;
}
void ClusterShapeAlgo::Create_Map ( const EcalRecHitCollection hits,
const CaloSubdetectorTopology topology 
) [private]

Definition at line 111 of file ClusterShapeAlgo.cc.

References eMaxId_, edm::SortedCollection< T, SORT >::end(), CaloRecHit::energy(), energyMap_, edm::SortedCollection< T, SORT >::find(), CaloNavigator< T >::home(), EcalRecHit::id(), CaloNavigator< T >::offsetBy(), x, and detailsBasic3DVector::y.

Referenced by Calculate().

{
  EcalRecHit tempEcalRecHit;
  CaloNavigator<DetId> posCurrent = CaloNavigator<DetId>(eMaxId_,topology );

  for(int x = 0; x < 5; x++)
    for(int y = 0; y < 5; y++)
    {
      posCurrent.home();
      posCurrent.offsetBy(-2+x,-2+y);

      if((*posCurrent != DetId(0)) && (hits->find(*posCurrent) != hits->end()))
      {
                                EcalRecHitCollection::const_iterator itt = hits->find(*posCurrent);
                                tempEcalRecHit = *itt;
                                energyMap_[y][x] = std::make_pair(tempEcalRecHit.id(),tempEcalRecHit.energy()); 
      }
      else
                                energyMap_[y][x] = std::make_pair(DetId(0), 0);  
    }
}
double ClusterShapeAlgo::f00 ( double  r) [private]

Definition at line 512 of file ClusterShapeAlgo.cc.

Referenced by Calculate_Polynomials().

{ return 1; }
double ClusterShapeAlgo::f11 ( double  r) [private]

Definition at line 514 of file ClusterShapeAlgo.cc.

References csvReporter::r.

Referenced by Calculate_Polynomials().

{ return r; }
double ClusterShapeAlgo::f20 ( double  r) [private]

Definition at line 516 of file ClusterShapeAlgo.cc.

Referenced by Calculate_Polynomials().

{ return 2.0*r*r-1.0; }
double ClusterShapeAlgo::f22 ( double  r) [private]

Definition at line 518 of file ClusterShapeAlgo.cc.

References csvReporter::r.

Referenced by Calculate_Polynomials().

{ return r*r; }
double ClusterShapeAlgo::f31 ( double  r) [private]

Definition at line 520 of file ClusterShapeAlgo.cc.

References csvReporter::r.

Referenced by Calculate_Polynomials().

{ return 3.0*r*r*r - 2.0*r; }
double ClusterShapeAlgo::f33 ( double  r) [private]

Definition at line 522 of file ClusterShapeAlgo.cc.

References csvReporter::r.

Referenced by Calculate_Polynomials().

{ return r*r*r; }
double ClusterShapeAlgo::f40 ( double  r) [private]

Definition at line 524 of file ClusterShapeAlgo.cc.

Referenced by Calculate_Polynomials().

{ return 6.0*r*r*r*r-6.0*r*r+1.0; }
double ClusterShapeAlgo::f42 ( double  r) [private]

Definition at line 526 of file ClusterShapeAlgo.cc.

References csvReporter::r.

Referenced by Calculate_Polynomials().

{ return 4.0*r*r*r*r-3.0*r*r; }
double ClusterShapeAlgo::f44 ( double  r) [private]

Definition at line 528 of file ClusterShapeAlgo.cc.

References csvReporter::r.

Referenced by Calculate_Polynomials().

{ return r*r*r*r; }
double ClusterShapeAlgo::f51 ( double  r) [private]

Definition at line 530 of file ClusterShapeAlgo.cc.

References funct::pow(), and csvReporter::r.

Referenced by Calculate_Polynomials().

{ return 10.0*pow(r,5)-12.0*pow(r,3)+3.0*r; }
double ClusterShapeAlgo::f53 ( double  r) [private]

Definition at line 532 of file ClusterShapeAlgo.cc.

References funct::pow().

Referenced by Calculate_Polynomials().

{ return 5.0*pow(r,5) - 4.0*pow(r,3); }
double ClusterShapeAlgo::f55 ( double  r) [private]

Definition at line 534 of file ClusterShapeAlgo.cc.

References funct::pow().

Referenced by Calculate_Polynomials().

{ return pow(r,5); }
double ClusterShapeAlgo::factorial ( int  n) const [private]

Definition at line 677 of file ClusterShapeAlgo.cc.

References i, and n.

Referenced by calc_AbsZernikeMoment().

                                              {
  double res=1.0;
  for(int i=2; i<=n; i++) res*=(double) i;
  return res;
}
double ClusterShapeAlgo::fast_AbsZernikeMoment ( const reco::BasicCluster passedCluster,
int  n,
int  m,
double  R0 
) [private]

Definition at line 536 of file ClusterShapeAlgo.cc.

References Calculate_Polynomials(), funct::cos(), energyDistribution_, fcn_, i, getHLTprescales::index, phi, csvReporter::r, query::result, funct::sin(), and mathSSE::sqrt().

Referenced by absZernikeMoment().

                                                                        {
  double r,ph,e,Re=0,Im=0,result;
  double TotalEnergy = passedCluster.energy();
  int index = (n/2)*(n/2)+(n/2)+m;
  int clusterSize=energyDistribution_.size();
  if(clusterSize<3) return 0.0;

  for (int i=0; i<clusterSize; i++)
    { 
      r = energyDistribution_[i].r / R0;
      if (r<1) {
        fcn_.clear();
        Calculate_Polynomials(r);
        ph = (energyDistribution_[i]).phi;
        e = energyDistribution_[i].deposited_energy;
        Re = Re + e/TotalEnergy * fcn_[index] * cos( (double) m * ph);
        Im = Im - e/TotalEnergy * fcn_[index] * sin( (double) m * ph);
      }
    }
  result = sqrt(Re*Re+Im*Im);

  return result;
}

Member Data Documentation

double ClusterShapeAlgo::A20_ [private]

Definition at line 103 of file ClusterShapeAlgo.h.

Referenced by Calculate(), and Calculate_ComplexZernikeMoments().

double ClusterShapeAlgo::A42_ [private]

Definition at line 103 of file ClusterShapeAlgo.h.

Referenced by Calculate(), and Calculate_ComplexZernikeMoments().

double ClusterShapeAlgo::covEtaEta_ [private]

Definition at line 96 of file ClusterShapeAlgo.h.

Referenced by Calculate(), and Calculate_Covariances().

double ClusterShapeAlgo::covEtaPhi_ [private]

Definition at line 96 of file ClusterShapeAlgo.h.

Referenced by Calculate(), and Calculate_Covariances().

double ClusterShapeAlgo::covPhiPhi_ [private]

Definition at line 96 of file ClusterShapeAlgo.h.

Referenced by Calculate(), and Calculate_Covariances().

double ClusterShapeAlgo::e2nd_ [private]

Definition at line 97 of file ClusterShapeAlgo.h.

Referenced by Calculate(), and Calculate_2ndEnergy().

Definition at line 106 of file ClusterShapeAlgo.h.

Referenced by Calculate(), and Calculate_2ndEnergy().

double ClusterShapeAlgo::e2x2_ [private]

Definition at line 97 of file ClusterShapeAlgo.h.

Referenced by Calculate(), and Calculate_e2x2().

Definition at line 94 of file ClusterShapeAlgo.h.

Referenced by Calculate_e2x2(), and Calculate_e4x4().

Definition at line 94 of file ClusterShapeAlgo.h.

Referenced by Calculate_e2x2(), and Calculate_e4x4().

Definition at line 98 of file ClusterShapeAlgo.h.

Referenced by Calculate(), and Calculate_e2x5Bottom().

double ClusterShapeAlgo::e2x5Left_ [private]

Definition at line 98 of file ClusterShapeAlgo.h.

Referenced by Calculate(), and Calculate_e2x5Left().

double ClusterShapeAlgo::e2x5Right_ [private]

Definition at line 98 of file ClusterShapeAlgo.h.

Referenced by Calculate(), and Calculate_e2x5Right().

double ClusterShapeAlgo::e2x5Top_ [private]

Definition at line 98 of file ClusterShapeAlgo.h.

Referenced by Calculate(), and Calculate_e2x5Top().

double ClusterShapeAlgo::e3x2_ [private]

Definition at line 97 of file ClusterShapeAlgo.h.

Referenced by Calculate(), and Calculate_e3x2().

double ClusterShapeAlgo::e3x2Ratio_ [private]

Definition at line 99 of file ClusterShapeAlgo.h.

Referenced by Calculate(), and Calculate_e3x2().

double ClusterShapeAlgo::e3x3_ [private]

Definition at line 97 of file ClusterShapeAlgo.h.

Referenced by Calculate(), and Calculate_e3x3().

double ClusterShapeAlgo::e4x4_ [private]

Definition at line 97 of file ClusterShapeAlgo.h.

Referenced by Calculate(), and Calculate_e4x4().

double ClusterShapeAlgo::e5x5_ [private]

Definition at line 97 of file ClusterShapeAlgo.h.

Referenced by Calculate(), Calculate_Covariances(), and Calculate_e5x5().

double ClusterShapeAlgo::eMax_ [private]

Definition at line 97 of file ClusterShapeAlgo.h.

Referenced by Calculate(), and Calculate_TopEnergy().

std::vector<double> ClusterShapeAlgo::energyBasketFractionEta_ [private]

Definition at line 104 of file ClusterShapeAlgo.h.

Referenced by Calculate(), and Calculate_BarrelBasketEnergyFraction().

std::vector<double> ClusterShapeAlgo::energyBasketFractionPhi_ [private]

Definition at line 105 of file ClusterShapeAlgo.h.

Referenced by Calculate(), and Calculate_BarrelBasketEnergyFraction().

std::pair<DetId, double> ClusterShapeAlgo::energyMap_[5][5] [private]
double ClusterShapeAlgo::etaLat_ [private]

Definition at line 101 of file ClusterShapeAlgo.h.

Referenced by Calculate(), and Calculate_lat().

std::vector<double> ClusterShapeAlgo::fcn_ [private]

Definition at line 108 of file ClusterShapeAlgo.h.

Referenced by Calculate_Polynomials(), and fast_AbsZernikeMoment().

double ClusterShapeAlgo::lat_ [private]

Definition at line 100 of file ClusterShapeAlgo.h.

Referenced by Calculate(), and Calculate_lat().

Definition at line 91 of file ClusterShapeAlgo.h.

Referenced by Calculate_Covariances(), and Calculate_EnergyDepTopology().

double ClusterShapeAlgo::phiLat_ [private]

Definition at line 102 of file ClusterShapeAlgo.h.

Referenced by Calculate(), and Calculate_lat().