#include <PFClusterShapeAlgo.h>

Public Member Functions
reco::ClusterShapeCollection *	makeClusterShapes (edm::Handle< reco::PFClusterCollection > clusterHandle, edm::Handle< reco::PFRecHitCollection > rechitHandle, const CaloSubdetectorGeometry barrelGeo_p, const CaloSubdetectorTopology barrelTop_p, const CaloSubdetectorGeometry endcapGeo_p, const CaloSubdetectorTopology endcapTop_p)
	PFClusterShapeAlgo (bool useFractions, double w0)
	~PFClusterShapeAlgo ()
Private Types
enum	Direction { N, NE, E, SE, S, SW, W, NW }
enum	GeomTopoIndex { BARREL = 0, ENDCAP = 1 }
typedef std::map< DetId, RecHitWithFraction >	RecHitMap
Private Member Functions
double	addMapEnergies (int etaIndexLow, int etaIndexHigh, int phiIndexLow, int phiIndexHigh)
void	covariances ()
void	fill5x5Map ()
void	find_e2x2 ()
void	find_e2x5Bottom ()
void	find_e2x5Left ()
void	find_e2x5Right ()
void	find_e2x5Top ()
void	find_e3x2 ()
void	find_e3x3 ()
void	find_e4x4 ()
void	find_e5x5 ()
void	find_eMax_e2nd ()
int	findPFRHIndexFromDetId (unsigned int id)
const reco::PFRecHitFraction *	getFractionFromDetId (const DetId &id)
reco::ClusterShape	makeClusterShape ()
Private Attributes
double	covEtaEta_
double	covEtaPhi_
double	covPhiPhi_
reco::PFClusterRef	currentCluster_p
unsigned int	currentClusterIndex_
edm::Handle < reco::PFRecHitCollection >	currentRecHit_v_p
double	e2nd_
DetId	e2ndId_
double	e2x2_
double	e2x5Bottom_
double	e2x5Left_
double	e2x5Right_
double	e2x5Top_
double	e3x2_
double	e3x2Ratio_
double	e3x3_
double	e4x4_
double	e5x5_
double	eMax_
Direction	eMaxDir
DetId	eMaxId_
unsigned int	geomIndex
std::vector< const CaloSubdetectorGeometry * >	geomVector
RecHitWithFraction	map5x5 [5][5]
math::XYZVector	meanPosition_
unsigned int	topoIndex
std::vector< const CaloSubdetectorTopology * >	topoVector
double	totalE_
bool	useFractions_
double	w0_

Detailed Description

Definition at line 37 of file PFClusterShapeAlgo.h.

Member Typedef Documentation

typedef std::map<DetId, RecHitWithFraction> PFClusterShapeAlgo::RecHitMap [private]

Definition at line 39 of file PFClusterShapeAlgo.h.

Member Enumeration Documentation

enum PFClusterShapeAlgo::Direction [private]

Enumerator:

N
NE
E
SE
S
SW
W
NW

Definition at line 41 of file PFClusterShapeAlgo.h.

{ N, NE, E, SE, S, SW, W, NW };

enum PFClusterShapeAlgo::GeomTopoIndex [private]

Enumerator:

BARREL
ENDCAP

Definition at line 42 of file PFClusterShapeAlgo.h.

{ BARREL = 0, ENDCAP = 1 };

Constructor & Destructor Documentation

PFClusterShapeAlgo::PFClusterShapeAlgo	(	bool	useFractions,
		double	w0
	)		`[explicit]`

Definition at line 3 of file PFClusterShapeAlgo.cc.

References useFractions_, and w0_.

{
  useFractions_ = useFractions;
  w0_ = w0;
}

PFClusterShapeAlgo::~PFClusterShapeAlgo ( )

Definition at line 9 of file PFClusterShapeAlgo.cc.

{
}

Member Function Documentation

double PFClusterShapeAlgo::addMapEnergies	(	int	etaIndexLow,
		int	etaIndexHigh,
		int	phiIndexLow,
		int	phiIndexHigh
	)		`[private]`

Definition at line 208 of file PFClusterShapeAlgo.cc.

References RecHitWithFraction::energy, relval_parameters_module::energy, i, j, and map5x5.

Referenced by find_e2x2(), find_e2x5Bottom(), find_e2x5Left(), find_e2x5Right(), find_e2x5Top(), find_e3x2(), find_e3x3(), find_e4x4(), and find_e5x5().

{
  const int etaLow  = etaIndexLow  + 2;
  const int etaHigh = etaIndexHigh + 2;
  const int phiLow  = phiIndexLow  + 2;
  const int phiHigh = phiIndexHigh + 2;

  double energy = 0;

  for (int i = etaLow; i <= etaHigh; ++i)
    {
      for (int j = phiLow; j <= phiHigh; ++j)
        {
          energy += map5x5[i][j].energy;
        }
    }
  return energy;
}

void PFClusterShapeAlgo::covariances ( ) [private]

Definition at line 299 of file PFClusterShapeAlgo.cc.

References covEtaEta_, covEtaPhi_, covPhiPhi_, RecoTauValidation_cfi::denominator, dPhi(), relval_parameters_module::energy, i, j, create_public_lumi_plots::log, map5x5, max(), meanPosition_, Geom::pi(), position, totalE_, Geom::twoPi(), w(), and w0_.

Referenced by makeClusterShape().

{
  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)
        {
          const math::XYZVector & crystalPosition(map5x5[i][j].position);
          
          double dPhi = crystalPosition.phi() - meanPosition_.phi();
          if (dPhi > + Geom::pi()) { dPhi = Geom::twoPi() - dPhi; }
          if (dPhi < - Geom::pi()) { dPhi = Geom::twoPi() - dPhi; }

          const double dEta = crystalPosition.eta() - meanPosition_.eta();
          
          const double w = std::max(0.0, w0_ + log(map5x5[i][j].energy / totalE_));
          
          denominator += w;
          numeratorEtaEta += w * dEta * dEta;
          numeratorEtaPhi += w * dEta * dPhi;
          numeratorPhiPhi += w * dPhi * dPhi;
        }
    }

  covEtaEta_ = numeratorEtaEta / denominator;
  covEtaPhi_ = numeratorEtaPhi / denominator;
  covPhiPhi_ = numeratorPhiPhi / denominator;
}

void PFClusterShapeAlgo::fill5x5Map ( ) [private]

Definition at line 155 of file PFClusterShapeAlgo.cc.

References currentRecHit_v_p, RecHitWithFraction::detId, eMaxId_, RecHitWithFraction::energy, findPFRHIndexFromDetId(), reco::PFRecHitFraction::fraction(), getFractionFromDetId(), CaloNavigator< T >::home(), i, getHLTprescales::index, j, map5x5, meanPosition_, CaloNavigator< T >::offsetBy(), RecHitWithFraction::position, position, topoIndex, topoVector, totalE_, and useFractions_.

Referenced by makeClusterShape().

{
  // first get a navigator to the central element
  CaloNavigator<DetId> position = CaloNavigator<DetId>(eMaxId_, topoVector[topoIndex]);

  meanPosition_ = math::XYZVector(0.0, 0.0, 0.0);
  totalE_ = 0;

  for (int i = 0; i < 5; ++i)
    {
      for (int j = 0; j < 5; ++j)
        {
          position.home();
          position.offsetBy(i - 2, j - 2);

          RecHitWithFraction newEntry;
          newEntry.detId = DetId(0);
          newEntry.energy = 0.0;
          newEntry.position = math::XYZVector(0, 0, 0);

          if (*position != DetId(0)) // if this is a valid detId...
            {
              // ...find the corresponding PFRecHit index
              const int index = findPFRHIndexFromDetId((*position).rawId());

              if (index >= 0) // if a PFRecHit exists for this detId
                {
                  double fraction = 1.0;
                  if (useFractions_) // if the algorithm should use fractions
                    { 
                      fraction = 0.0;
                      const reco::PFRecHitFraction * fraction_p = getFractionFromDetId(*position);
                      if (fraction_p) { fraction = fraction_p->fraction(); }
                    }

                  const reco::PFRecHitRef rhRef(currentRecHit_v_p, index);
                  const math::XYZVector crystalPosition(rhRef->position());
                  const double energyFraction =  rhRef->energy() * fraction;

                  newEntry.detId = *position;
                  newEntry.energy = energyFraction;
                  newEntry.position = crystalPosition;

                  meanPosition_ = meanPosition_ + crystalPosition * energyFraction; 
                  totalE_ += energyFraction;
                }
            }
          map5x5[i][j] = newEntry;
        }
    }
  meanPosition_ /= totalE_;
}

void PFClusterShapeAlgo::find_e2x2 ( ) [private]

Definition at line 242 of file PFClusterShapeAlgo.cc.

References addMapEnergies(), e2x2_, eMaxDir, NE, NW, SE, and SW.

Referenced by makeClusterShape().

{
  std::map<double, Direction> directionMap;

  directionMap[addMapEnergies(-1, +0, -1, +0)] = SE;
  directionMap[addMapEnergies(-1, +0, +0, +1)] = NE;
  directionMap[addMapEnergies(+0, +1, -1, +0)] = SW;
  directionMap[addMapEnergies(+0, +1, +0, +1)] = NW;

  const std::map<double, Direction>::reverse_iterator eMaxDir_it = directionMap.rbegin();

  eMaxDir = eMaxDir_it->second;

  e2x2_ = eMaxDir_it->first;
}

void PFClusterShapeAlgo::find_e2x5Bottom ( ) [private]

Definition at line 232 of file PFClusterShapeAlgo.cc.

References addMapEnergies(), and e2x5Bottom_.

Referenced by makeClusterShape().

{ e2x5Bottom_ = addMapEnergies(+1, +2, -2, +2); }

void PFClusterShapeAlgo::find_e2x5Left ( ) [private]

Definition at line 230 of file PFClusterShapeAlgo.cc.

References addMapEnergies(), and e2x5Left_.

Referenced by makeClusterShape().

{ e2x5Left_   = addMapEnergies(-2, +2, -2, -1); }

void PFClusterShapeAlgo::find_e2x5Right ( ) [private]

Definition at line 229 of file PFClusterShapeAlgo.cc.

References addMapEnergies(), and e2x5Right_.

Referenced by makeClusterShape().

{ e2x5Right_  = addMapEnergies(-2, +2, +1, +2); }

void PFClusterShapeAlgo::find_e2x5Top ( ) [private]

Definition at line 231 of file PFClusterShapeAlgo.cc.

References addMapEnergies(), and e2x5Top_.

Referenced by makeClusterShape().

{ e2x5Top_    = addMapEnergies(-2, -1, -2, +2); }

void PFClusterShapeAlgo::find_e3x2 ( ) [private]

Definition at line 258 of file PFClusterShapeAlgo.cc.

References addMapEnergies(), RecoTauValidation_cfi::denominator, dir, E, e3x2_, e3x2Ratio_, RecHitWithFraction::energy, map5x5, N, S, and W.

Referenced by makeClusterShape().

{
  // Find the direction of the highest energy neighbour
  std::map<double, Direction> directionMap;
  directionMap[map5x5[2][3].energy] = N; 
  directionMap[map5x5[2][1].energy] = S;
  directionMap[map5x5[1][2].energy] = E;
  directionMap[map5x5[3][2].energy] = W;
  // Maps are sorted in ascending order - get the last element
  const Direction dir = directionMap.rbegin()->second;

  if (dir == N) 
    {
      e3x2_ = addMapEnergies(-1, +1, +0, +1);
      const double numerator   = map5x5[3][2].energy + map5x5[1][2].energy;
      const double denominator = map5x5[1][3].energy + map5x5[3][3].energy + 0.5;
      e3x2Ratio_ = numerator / denominator;
    }
  else if (dir == S)
    {
      e3x2_ = addMapEnergies(-1, +1, -1, +0);
      const double numerator   = map5x5[3][2].energy + map5x5[1][2].energy;
      const double denominator = map5x5[1][1].energy + map5x5[3][1].energy + 0.5;
      e3x2Ratio_ = numerator / denominator;
    }
  else if (dir == W)
    {
      e3x2_ = addMapEnergies(+0, +1, -1, +1);
      const double numerator   = map5x5[2][3].energy + map5x5[2][1].energy;
      const double denominator = map5x5[3][3].energy + map5x5[3][1].energy + 0.5;
      e3x2Ratio_ = numerator / denominator;
    }
  else if (dir == E)
    {
      e3x2_ = addMapEnergies(-1, +0, -1, +1);
      const double numerator   = map5x5[2][3].energy + map5x5[2][1].energy;
      const double denominator = map5x5[1][1].energy + map5x5[1][3].energy + 0.5;
      e3x2Ratio_ = numerator / denominator;
    }
}

void PFClusterShapeAlgo::find_e3x3 ( ) [private]

Definition at line 227 of file PFClusterShapeAlgo.cc.

References addMapEnergies(), and e3x3_.

Referenced by makeClusterShape().

{ e3x3_       = addMapEnergies(-1, +1, -1, +1); }

void PFClusterShapeAlgo::find_e4x4 ( ) [private]

Definition at line 234 of file PFClusterShapeAlgo.cc.

References addMapEnergies(), e4x4_, eMaxDir, NE, NW, SE, and SW.

Referenced by makeClusterShape().

{
  if (eMaxDir == SE) { e4x4_ = addMapEnergies(-2, +1, -2, +1); return; }
  if (eMaxDir == NE) { e4x4_ = addMapEnergies(-2, +1, -1, +2); return; }
  if (eMaxDir == SW) { e4x4_ = addMapEnergies(-1, +2, -2, +1); return; }
  if (eMaxDir == NW) { e4x4_ = addMapEnergies(-1, +2, -1, +2); return; }
}

void PFClusterShapeAlgo::find_e5x5 ( ) [private]

Definition at line 228 of file PFClusterShapeAlgo.cc.

References addMapEnergies(), and e5x5_.

Referenced by makeClusterShape().

{ e5x5_       = addMapEnergies(-2, +2, -2, +2); }

void PFClusterShapeAlgo::find_eMax_e2nd ( ) [private]

Definition at line 96 of file PFClusterShapeAlgo.cc.

References currentCluster_p, e2nd_, e2ndId_, eMax_, and eMaxId_.

Referenced by makeClusterShape().

{
  std::map<double, DetId> energyMap;

  // First get the RecHitFractions:
  const std::vector<reco::PFRecHitFraction> & fraction_v = currentCluster_p->recHitFractions();
  // For every one of them...
  for (std::vector<reco::PFRecHitFraction>::const_iterator it = fraction_v.begin(); it != fraction_v.end(); ++it)
    {
      // ...find the corresponding rechit:
      // const reco::PFRecHit & rechit = (*currentRecHit_v_p)[it->recHitIndex()];
      const reco::PFRecHitRef rechit = it->recHitRef();
      // ...and DetId:
      const DetId rechitDetId = DetId(rechit->detId());
      // Make the new Pair and put it in the map:
      energyMap[rechit->energy()] = rechitDetId;
    }
  // maps are sorted in ascending order so get the last two elements:
  std::map<double, DetId>::reverse_iterator it = energyMap.rbegin();
  eMax_   = it->first;
  eMaxId_ = it->second;
  it++;
  e2nd_   = it->first;
  e2ndId_ = it->second;
}

int PFClusterShapeAlgo::findPFRHIndexFromDetId ( unsigned int id ) [private]

Definition at line 122 of file PFClusterShapeAlgo.cc.

References currentRecHit_v_p, errorMatrix2Lands_multiChannel::id, getHLTprescales::index, and gen::k.

Referenced by fill5x5Map().

{
  int index = -1; // need some negative number
  for (unsigned int k = 0; k < currentRecHit_v_p->size(); ++k)
    {
      if ((*currentRecHit_v_p)[k].detId() == id)
        {
          index = static_cast<int>(k);
          break;
        }
    }
  return index;
}

const reco::PFRecHitFraction * PFClusterShapeAlgo::getFractionFromDetId ( const DetId & id ) [private]

Definition at line 137 of file PFClusterShapeAlgo.cc.

References currentCluster_p.

Referenced by fill5x5Map().

{
  const std::vector< reco::PFRecHitFraction > & fraction_v = currentCluster_p->recHitFractions();
  for (std::vector<reco::PFRecHitFraction>::const_iterator it = fraction_v.begin(); it != fraction_v.end(); ++it)
    {
      //const unsigned int rhIndex = it->recHitIndex();
      //reco::PFRecHitRef rh_p(currentRecHit_v_p, rhIndex);
      const reco::PFRecHitRef rh_p = it->recHitRef();
      const DetId rhDetId = DetId(rh_p->detId());
      if (rhDetId == id) 
        { 
          return &(*it); 
        }
    }
  return 0;
}

reco::ClusterShape PFClusterShapeAlgo::makeClusterShape ( ) [private]

Definition at line 60 of file PFClusterShapeAlgo.cc.

References covariances(), covEtaEta_, covEtaPhi_, covPhiPhi_, e2nd_, e2ndId_, e2x2_, e2x5Bottom_, e2x5Left_, e2x5Right_, e2x5Top_, e3x2_, e3x2Ratio_, e3x3_, e4x4_, e5x5_, eMax_, eMaxId_, fill5x5Map(), find_e2x2(), find_e2x5Bottom(), find_e2x5Left(), find_e2x5Right(), find_e2x5Top(), find_e3x2(), find_e3x3(), find_e4x4(), find_e5x5(), and find_eMax_e2nd().

Referenced by makeClusterShapes().

{
  find_eMax_e2nd();
  fill5x5Map();
  
  find_e2x2();
  find_e3x2();
  find_e3x3();
  find_e4x4();
  find_e5x5();
  
  find_e2x5Right();
  find_e2x5Left();
  find_e2x5Top();
  find_e2x5Bottom();
  
  covariances();
  
  double dummyLAT = 0;
  double dummyEtaLAT = 0;
  double dummyPhiLAT = 0;
  double dummyA20 = 0;
  double dummyA42 = 0;

  std::vector<double> dummyEnergyBasketFractionEta_v;
  std::vector<double> dummyEnergyBasketFractionPhi_v;

  return reco::ClusterShape(covEtaEta_, covEtaPhi_, covPhiPhi_, 
                            eMax_, eMaxId_, e2nd_, e2ndId_,
                            e2x2_, e3x2_, e3x3_, e4x4_, e5x5_, 
                            e2x5Right_, e2x5Left_, e2x5Top_, e2x5Bottom_, e3x2Ratio_,
                            dummyLAT, dummyEtaLAT, dummyPhiLAT, dummyA20, dummyA42,
                            dummyEnergyBasketFractionEta_v, dummyEnergyBasketFractionPhi_v);
}

reco::ClusterShapeCollection * PFClusterShapeAlgo::makeClusterShapes	(	edm::Handle< reco::PFClusterCollection >	clusterHandle,
		edm::Handle< reco::PFRecHitCollection >	rechitHandle,
		const CaloSubdetectorGeometry *	barrelGeo_p,
		const CaloSubdetectorTopology *	barrelTop_p,
		const CaloSubdetectorGeometry *	endcapGeo_p,
		const CaloSubdetectorTopology *	endcapTop_p
	)

Definition at line 14 of file PFClusterShapeAlgo.cc.

References BARREL, currentCluster_p, currentClusterIndex_, currentRecHit_v_p, ENDCAP, geomIndex, geomVector, i, makeClusterShape(), topoIndex, and topoVector.

{
  static const float etaEndOfBarrel = 1.497;

  topoVector.push_back(the_barrelTop_p);
  topoVector.push_back(the_endcapTop_p);
  geomVector.push_back(the_barrelGeo_p);
  geomVector.push_back(the_endcapGeo_p);

  reco::ClusterShapeCollection * shape_v_p = new reco::ClusterShapeCollection();

  currentRecHit_v_p = rechitHandle;

  for (unsigned int i = 0; i < clusterHandle->size(); ++i)
    {
      // Make each cluster the "current" cluster
      currentCluster_p = reco::PFClusterRef(clusterHandle, i);
      currentClusterIndex_ = i;

      // Find the right topology to use with this cluster
      topoIndex = BARREL;
      geomIndex = BARREL;
      const math::XYZVector currentClusterPos(currentCluster_p->position());
      if (fabs(currentClusterPos.eta()) > etaEndOfBarrel)
        {
          topoIndex = ENDCAP;
          geomIndex = ENDCAP;
        }

      // Create the clustershape and push it into the vector
      shape_v_p->push_back(makeClusterShape());
    }

  topoVector.clear();
  topoVector.clear();
  geomVector.clear();
  geomVector.clear();

  return shape_v_p;
}