#include <PositionCalc.h>

Public Types
typedef std::vector< std::pair < DetId, double > >	HitsAndEnergies

typedef std::vector< std::pair < DetId, float > >	HitsAndFractions

Public Member Functions
template<typename HitType >
math::XYZPoint	Calculate_Location (const HitsAndFractions &iDetIds, const edm::SortedCollection< HitType > iRecHits, const CaloSubdetectorGeometry iSubGeom, const CaloSubdetectorGeometry *iESGeom=0)

const PositionCalc &	operator= (const PositionCalc &rhs)

	PositionCalc (const edm::ParameterSet &par)

	PositionCalc ()

Private Attributes
const CaloSubdetectorGeometry *	m_esGeom

bool	m_esMinus

bool	m_esPlus

bool	param_LogWeighted_

double	param_T0_barl_

double	param_T0_endc_

double	param_T0_endcPresh_

double	param_W0_

double	param_X0_

Detailed Description

Definition at line 29 of file PositionCalc.h.

Member Typedef Documentation

typedef std::vector< std::pair<DetId,double> > PositionCalc::HitsAndEnergies

Definition at line 33 of file PositionCalc.h.

typedef std::vector< std::pair<DetId,float> > PositionCalc::HitsAndFractions

Definition at line 32 of file PositionCalc.h.

Constructor & Destructor Documentation

PositionCalc::PositionCalc ( const edm::ParameterSet & par )

Definition at line 4 of file PositionCalc.cc.

                                                      :
   param_LogWeighted_   ( par.getParameter<bool>("LogWeighted")) ,
   param_T0_barl_       ( par.getParameter<double>("T0_barl")) , 
   param_T0_endc_       ( par.getParameter<double>("T0_endc")) , 
   param_T0_endcPresh_  ( par.getParameter<double>("T0_endcPresh")) , 
   param_W0_            ( par.getParameter<double>("W0")) ,
   param_X0_            ( par.getParameter<double>("X0")) ,
   m_esGeom             ( 0 ) ,
   m_esPlus             ( false ) ,
   m_esMinus            ( false )
 {
 }

PositionCalc::PositionCalc ( )

inline

Definition at line 38 of file PositionCalc.h.

38 { };

Member Function Documentation

template<typename HitType >

math::XYZPoint PositionCalc::Calculate_Location	(	const HitsAndFractions &	iDetIds,
		const edm::SortedCollection< HitType > *	iRecHits,
		const CaloSubdetectorGeometry *	iSubGeom,
		const CaloSubdetectorGeometry *	iESGeom = `0`
	)

Definition at line 68 of file PositionCalc.h.

References EcalBarrel, edm::SortedCollection< T, SORT >::end(), relval_parameters_module::energy, edm::hlt::Exception, edm::SortedCollection< T, SORT >::find(), cropTnPTrees::frac, CaloSubdetectorGeometry::getGeometry(), CaloSubdetectorGeometry::getValidDetIds(), j, log, LogDebug, m_esGeom, m_esMinus, m_esPlus, bookConverter::max, n, DetId::null(), param_LogWeighted_, param_T0_barl_, param_T0_endc_, param_T0_endcPresh_, param_W0_, param_X0_, edm::SortedCollection< T, SORT >::size(), DetId::subdetId(), histoStyle::weight, detailsBasic3DVector::z, and zw().

Referenced by HLTAlCaMonPi0::analyze(), DQMSourcePi0::analyze(), HLTAlCa::analyze(), HybridClusterAlgo::e2Et(), HybridClusterAlgo::et25(), PhotonProducer::fillPhotonCollection(), Multi5x5ClusterAlgo::mainSearch(), HybridClusterAlgo::mainSearch(), IslandClusterAlgo::makeCluster(), CosmicClusterAlgo::makeCluster(), Multi5x5ClusterAlgo::makeCluster(), and EgammaHLTNxNClusterProducer::makeNxNClusters().

                                                              {
   typedef edm::SortedCollection<HitType> HitTypeCollection;
   math::XYZPoint returnValue ( 0, 0, 0 ) ;
   
   // Throw an error if the cluster was not initialized properly
   
   if( 0 == iRecHits || 0 == iSubGeom ) {
     throw cms::Exception("PositionCalc")
       << "Calculate_Location() called uninitialized or wrong initialization.";
   }
   
   if( 0 != iDetIds.size()   &&
       0 != iRecHits->size()     ) {
     
     HitsAndEnergies detIds; 
     detIds.reserve( iDetIds.size() ) ;
     
     double eTot  ( 0 ) ;
     double eMax  ( 0 ) ;
     DetId  maxId ;
     
     // Check that DetIds are nonzero
     typename HitTypeCollection::const_iterator endRecHits ( iRecHits->end() ) ;
     HitsAndFractions::const_iterator n, endDiDs( iDetIds.end() );
     for( n = iDetIds.begin(); n != endDiDs ; ++n ) {
       const DetId dId ( (*n).first ) ;
       const float frac( (*n).second) ; 
       if( !dId.null() ) {
     typename HitTypeCollection::const_iterator iHit ( iRecHits->find( dId ) ) ;
     if( iHit != endRecHits ) {                 
       const double energy ( iHit->energy() *frac ) ;       
       detIds.push_back( std::make_pair(dId,energy) );
       if( 0.0 < energy ) { // only save positive energies       
         if( eMax < energy ) {
           eMax  = energy ;
           maxId = dId    ;        
         }
         eTot += energy ;
       }
     }
       }
     }
     
     if( 0.0 >= eTot ) {
       LogDebug("ZeroClusterEnergy") << "cluster with 0 energy: " 
                     << eTot << " size: " << detIds.size() 
                     << " , returning (0,0,0)";
     } else {
       // first time or when es geom changes set flags
       if( 0 != iESGeom && m_esGeom != iESGeom ) {
     m_esGeom = iESGeom ;
     for( uint32_t ic ( 0 ) ;
          ( ic != m_esGeom->getValidDetIds().size() ) &&
            ( (!m_esPlus) || (!m_esMinus) ) ; ++ic ) {
       const double z ( m_esGeom->getGeometry( m_esGeom->getValidDetIds()[ic] )->getPosition().z() ) ;
       m_esPlus  = m_esPlus  || ( 0 < z ) ;
       m_esMinus = m_esMinus || ( 0 > z ) ;
     }
       }
       
       //Select the correct value of the T0 parameter depending on subdetector       
       const CaloCellGeometry* center_cell ( iSubGeom->getGeometry( maxId ) ) ;
       const double ctreta (center_cell->getPosition().eta());
       
       // for barrel, use barrel T0; 
       // for endcap: if preshower present && in preshower fiducial, 
       //             use preshower T0
       // else use endcap only T0       
       const double preshowerStartEta =  1.653;
       const int subdet = maxId.subdetId();
       const double T0 ( subdet == EcalBarrel ? param_T0_barl_ :
             ( ( ( preshowerStartEta < fabs( ctreta ) ) &&
                 ( ( ( 0 < ctreta ) && 
                 m_esPlus          ) ||
                   ( ( 0 > ctreta ) &&
                 m_esMinus         )   )    ) ?
               param_T0_endcPresh_ : param_T0_endc_ ) ) ;
       
       // Calculate shower depth
       const float maxDepth ( param_X0_ * ( T0 + log( eTot ) ) ) ;
       const float maxToFront ( center_cell->getPosition().mag() ) ; // to front face
       
       // Loop over hits and get weights
       double total_weight = 0;
       const double eTot_inv = 1.0/eTot;
       const double logETot_inv = ( param_LogWeighted_ ? log(eTot_inv) : 0 );
       
       double xw ( 0 ) ;
       double yw ( 0 ) ;
       double zw ( 0 ) ;
       
       
       HitsAndEnergies::const_iterator j, hAndE_end = detIds.end();
       for( j = detIds.begin() ; j != hAndE_end ; ++j ) {
     const DetId dId ( (*j).first )  ;
     const double e_j( (*j).second ) ;        
     
     double weight = 0;
     if ( param_LogWeighted_ ) {
       if ( e_j > 0.0 ) {
         weight = std::max( 0., param_W0_ + log(e_j) + logETot_inv );
       } else {
         weight = 0;
       }
     } else {
       weight = e_j*eTot_inv;
     }
     
     const CaloCellGeometry* cell ( iSubGeom->getGeometry( dId ) ) ;
     const float depth ( maxDepth + maxToFront - cell->getPosition().mag() ) ;
     
     const GlobalPoint pos (static_cast<const TruncatedPyramid*>( cell )->getPosition( depth ) );
     
     xw += weight*pos.x() ;
     yw += weight*pos.y() ;
     zw += weight*pos.z() ;
     
     total_weight += weight ;
       }
       returnValue = math::XYZPoint( xw/total_weight, 
                     yw/total_weight, 
                     zw/total_weight ) ;
     }
   }
   return returnValue ;
 }

const PositionCalc & PositionCalc::operator= ( const PositionCalc & rhs )

Definition at line 17 of file PositionCalc.cc.

References m_esGeom, m_esMinus, m_esPlus, param_LogWeighted_, param_T0_barl_, param_T0_endc_, param_T0_endcPresh_, param_W0_, and param_X0_.

 {
    param_LogWeighted_ = rhs.param_LogWeighted_;
    param_T0_barl_ = rhs.param_T0_barl_;
    param_T0_endc_ = rhs.param_T0_endc_;
    param_T0_endcPresh_ = rhs.param_T0_endcPresh_;
    param_W0_ = rhs.param_W0_;
    param_X0_ = rhs.param_X0_;
 
    m_esGeom = rhs.m_esGeom ;
    m_esPlus = rhs.m_esPlus ;
    m_esMinus = rhs.m_esMinus ;
    return *this;
 }