#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=nullptr)

const PositionCalc &	operator= (const PositionCalc &rhs)

	PositionCalc ()

	PositionCalc (const edm::ParameterSet &par)

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

◆ HitsAndEnergies

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

Definition at line 32 of file PositionCalc.h.

◆ HitsAndFractions

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

Definition at line 31 of file PositionCalc.h.

Constructor & Destructor Documentation

◆ PositionCalc() [1/2]

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(nullptr),
       m_esPlus(false),
       m_esMinus(false) {}

◆ PositionCalc() [2/2]

PositionCalc::PositionCalc ( )

inline

Definition at line 37 of file PositionCalc.h.

37 {};

Member Function Documentation

◆ Calculate_Location()

template<typename HitType >

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

Definition at line 65 of file PositionCalc.h.

                                                                                         {
   typedef edm::SortedCollection<HitType> HitTypeCollection;
   math::XYZPoint returnValue(0, 0, 0);
  
   // Throw an error if the cluster was not initialized properly
  
   if (nullptr == iRecHits || nullptr == iSubGeom) {
     throw cms::Exception("PositionCalc") << "Calculate_Location() called uninitialized or wrong initialization.";
   }
  
   if (!iDetIds.empty() && !iRecHits->empty()) {
     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 (nullptr != 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
       auto 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;
         }
  
         auto cell(iSubGeom->getGeometry(dId));
         const float depth(maxDepth + maxToFront - cell->getPosition().mag());
  
         const GlobalPoint pos(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;
 }

References LEDCalibrationChannels::depth, EcalBarrel, cosmicPhotonAnalyzer_cfi::eMax, edm::SortedCollection< T, SORT >::empty(), edm::SortedCollection< T, SORT >::end(), HCALHighEnergyHPDFilter_cfi::energy, Exception, edm::SortedCollection< T, SORT >::find(), cropTnPTrees::frac, CaloSubdetectorGeometry::getGeometry(), CaloSubdetectorGeometry::getValidDetIds(), dqmiolumiharvest::j, dqm-mbProfile::log, LogDebug, m_esGeom, m_esMinus, m_esPlus, SiStripPI::max, HLT_2018_cff::maxDepth, dqmiodumpmetadata::n, DetId::null(), param_LogWeighted_, param_T0_barl_, param_T0_endc_, param_T0_endcPresh_, param_W0_, param_X0_, DetId::subdetId(), mps_merge::weight, z, and zw().

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

◆ operator=()

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

Definition at line 15 of file PositionCalc.cc.

                                                                    {
   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;
 }