#include <PositionCalc.h>

Public Member Functions
math::XYZPoint	Calculate_Location (const std::vector< std::pair< DetId, float > > &iDetIds, const EcalRecHitCollection 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 25 of file PositionCalc.h.

Constructor & Destructor Documentation

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

Definition at line 11 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 32 of file PositionCalc.h.

32 { };

Member Function Documentation

math::XYZPoint PositionCalc::Calculate_Location	(	const std::vector< std::pair< DetId, float > > &	iDetIds,
		const EcalRecHitCollection *	iRecHits,
		const CaloSubdetectorGeometry *	iSubGeom,
		const CaloSubdetectorGeometry *	iESGeom = `0`
	)

Definition at line 40 of file PositionCalc.cc.

References EcalBarrel, jptDQMConfig_cff::eMax, edm::SortedCollection< T, SORT >::end(), relval_parameters_module::energy, edm::SortedCollection< T, SORT >::find(), CaloSubdetectorGeometry::getGeometry(), CaloSubdetectorGeometry::getValidDetIds(), j, create_public_lumi_plots::log, LogDebug, m_esGeom, m_esMinus, m_esPlus, max(), n, DetId::null(), param_LogWeighted_, param_T0_barl_, param_T0_endc_, param_T0_endcPresh_, param_W0_, param_X0_, pos, edm::SortedCollection< T, SORT >::size(), DetId::subdetId(), CommonMethods::weight(), and detailsBasic3DVector::z.

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

 {
    math::XYZPoint returnValue ( 0, 0, 0 ) ;
 
    // Throw an error if the cluster was not initialized properly
 
    if( 0 == iRecHits || 
        0 == iSubGeom    )
    {
       throw(std::runtime_error("\n\nPositionCalc::Calculate_Location called uninitialized or wrong initialization.\n\n"));
    }
 
    if( 0 != iDetIds.size()   &&
        0 != iRecHits->size()     )
    {
       typedef std::vector<DetId> DetIdVec ;
 
       DetIdVec detIds ;
       detIds.reserve( iDetIds.size() ) ;
 
       double eTot  ( 0 ) ;
       double eMax  ( 0 ) ;
       DetId  maxId ;
 
       // Check that DetIds are nonzero
       EcalRecHitCollection::const_iterator endRecHits ( iRecHits->end() ) ;
       for( std::vector< std::pair<DetId, float> >::const_iterator n ( iDetIds.begin() ) ; n != iDetIds.end() ; ++n ) 
       {
      const DetId dId ( (*n).first ) ;
      if( !dId.null() )
      {
         EcalRecHitCollection::const_iterator iHit ( iRecHits->find( dId ) ) ;
         if( iHit != endRecHits )
         {
            detIds.push_back( dId );
 
            const double energy ( iHit->energy() ) ;
 
            if( 0 < energy ) // only save positive energies
            {
           if( eMax < energy )
           {
              eMax  = energy ;
              maxId = dId    ;
           }
           eTot += energy ;
            }
         }
      }
       }
 
       if( 0 >= eTot )
       {
      LogDebug("ZeroClusterEnergy") << "cluster with 0 energy: " << eTot
                        << ", 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;
 
      double xw ( 0 ) ;
      double yw ( 0 ) ;
      double zw ( 0 ) ;
 
      for( DetIdVec::const_iterator j ( detIds.begin() ) ; j != detIds.end() ; ++j ) 
      {
         const DetId dId ( *j ) ;
         EcalRecHitCollection::const_iterator iR ( iRecHits->find( dId ) ) ;
         const double e_j ( iR->energy() ) ;
 
         double weight = 0;
             if ( param_LogWeighted_ ) {
                     if ( e_j > 0 ) {
                             weight = std::max( 0., param_W0_ + log(e_j/eTot) );
                     } else {
                             weight = 0;
                     }
             } else {
                     weight = e_j/eTot;
             }
 
         const CaloCellGeometry* cell ( iSubGeom->getGeometry( dId ) ) ;
         const float depth ( maxDepth + maxToFront - cell->getPosition().mag() ) ;
 
         const GlobalPoint pos (
            dynamic_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 24 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;
 }