This is EGM version of the ECAL position + depth correction calculation. More...

Inheritance diagram for ECAL2DPositionCalcWithDepthCorr:

Public Member Functions
void	calculateAndSetPosition (reco::PFCluster &, const HcalPFCuts *) override

void	calculateAndSetPositions (reco::PFClusterCollection &, const HcalPFCuts *) override

	ECAL2DPositionCalcWithDepthCorr (const edm::ParameterSet &conf, edm::ConsumesCollector &cc)

	ECAL2DPositionCalcWithDepthCorr (const ECAL2DPositionCalcWithDepthCorr &)=delete

ECAL2DPositionCalcWithDepthCorr &	operator= (const ECAL2DPositionCalcWithDepthCorr &)=delete

void	update (const edm::EventSetup &es) override

Public Member Functions inherited from PFCPositionCalculatorBase
const std::string &	name () const

PosCalc &	operator= (const PosCalc &)=delete

	PFCPositionCalculatorBase (const edm::ParameterSet &conf, edm::ConsumesCollector &cc)

	PFCPositionCalculatorBase (const PosCalc &)=delete

virtual	~PFCPositionCalculatorBase ()=default

Private Member Functions
void	calculateAndSetPositionActual (reco::PFCluster &) const

Private Attributes
const CaloSubdetectorGeometry *	_ebGeom

const CaloSubdetectorGeometry *	_eeGeom

const CaloSubdetectorGeometry *	_esGeom

bool	_esMinus

bool	_esPlus

edm::ESGetToken< CaloGeometry, CaloGeometryRecord >	_geomToken

const double	_minAllowedNorm

const double	_param_T0_EB

const double	_param_T0_EE

const double	_param_T0_ES

const double	_param_W0

const double	_param_X0

std::unique_ptr< CaloRecHitResolutionProvider >	_timeResolutionCalc

Additional Inherited Members
Protected Attributes inherited from PFCPositionCalculatorBase
const float	_minFractionInCalc

Detailed Description

This is EGM version of the ECAL position + depth correction calculation.

Definition at line 24 of file ECAL2DPositionCalcWithDepthCorr.cc.

Constructor & Destructor Documentation

◆ ECAL2DPositionCalcWithDepthCorr() [1/2]

ECAL2DPositionCalcWithDepthCorr::ECAL2DPositionCalcWithDepthCorr	(	const edm::ParameterSet &	conf,
		edm::ConsumesCollector &	cc
	)

inline

Definition at line 26 of file ECAL2DPositionCalcWithDepthCorr.cc.

References _timeResolutionCalc, edm::ParameterSet::exists(), and edm::ParameterSet::getParameterSet().

       : PFCPositionCalculatorBase(conf, cc),
         _param_T0_EB(conf.getParameter<double>("T0_EB")),
         _param_T0_EE(conf.getParameter<double>("T0_EE")),
         _param_T0_ES(conf.getParameter<double>("T0_ES")),
         _param_W0(conf.getParameter<double>("W0")),
         _param_X0(conf.getParameter<double>("X0")),
         _minAllowedNorm(conf.getParameter<double>("minAllowedNormalization")),
         _ebGeom(nullptr),
         _eeGeom(nullptr),
         _esGeom(nullptr),
         _esPlus(false),
         _esMinus(false),
         _geomToken(cc.esConsumes<edm::Transition::BeginRun>()) {
     _timeResolutionCalc.reset(nullptr);
     if (conf.exists("timeResolutionCalc")) {
       const edm::ParameterSet& timeResConf = conf.getParameterSet("timeResolutionCalc");
       _timeResolutionCalc = std::make_unique<CaloRecHitResolutionProvider>(timeResConf);
     }
   }

◆ ECAL2DPositionCalcWithDepthCorr() [2/2]

ECAL2DPositionCalcWithDepthCorr::ECAL2DPositionCalcWithDepthCorr ( const ECAL2DPositionCalcWithDepthCorr & )

delete

Member Function Documentation

◆ calculateAndSetPosition()

void ECAL2DPositionCalcWithDepthCorr::calculateAndSetPosition	(	reco::PFCluster &	cluster,
		const HcalPFCuts *
	)

overridevirtual

Implements PFCPositionCalculatorBase.

Definition at line 95 of file ECAL2DPositionCalcWithDepthCorr.cc.

References calculateAndSetPositionActual().

                                                                                                        {
   calculateAndSetPositionActual(cluster);
 }

◆ calculateAndSetPositionActual()

void ECAL2DPositionCalcWithDepthCorr::calculateAndSetPositionActual ( reco::PFCluster & cluster ) const

private

Definition at line 105 of file ECAL2DPositionCalcWithDepthCorr.cc.

References _ebGeom, _eeGeom, _esMinus, _esPlus, _minAllowedNorm, _param_T0_EB, _param_T0_EE, _param_T0_ES, _param_W0, _param_X0, _timeResolutionCalc, funct::abs(), reco::PFCluster::calculatePositionREP(), ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), hcalRecHitTable_cff::depth, PFLayer::ECAL_BARREL, PFLayer::ECAL_ENDCAP, reco::PFCluster::energy(), Exception, f, nano_mu_digi_cff::float, CaloSubdetectorGeometry::getGeometry(), edm::isFinite(), SiStripPI::max, HLT_2024v11_cff::maxDepth, PFLayer::NONE, reco::PFCluster::recHitFractions(), reco::CaloCluster::seed(), reco::CaloCluster::setEnergy(), reco::PFCluster::setLayer(), reco::CaloCluster::setPosition(), reco::PFCluster::setTime(), reco::PFCluster::setTimeError(), mathSSE::sqrt(), mps_merge::weight, x, y, and z.

Referenced by calculateAndSetPosition(), and calculateAndSetPositions().

                                                                                                 {
   constexpr double preshowerStartEta = 1.653;
   constexpr double preshowerEndEta = 2.6;
   if (!cluster.seed()) {
     throw cms::Exception("ClusterWithNoSeed") << " Found a cluster with no seed: " << cluster;
   }
   double cl_energy = 0;
   double cl_energy_float = 0;
   double cl_time = 0;
   double cl_timeweight = 0.0;
   double max_e = 0.0;
   double clusterT0 = 0.0;
   PFLayer::Layer max_e_layer = PFLayer::NONE;
   reco::PFRecHitRef refmax;
   // find the seed and max layer
   for (const reco::PFRecHitFraction& rhf : cluster.recHitFractions()) {
     const reco::PFRecHitRef& refhit = rhf.recHitRef();
     const double rh_fraction = rhf.fraction();
     const double rh_rawenergy = refhit->energy();
     const double rh_energy = rh_rawenergy * rh_fraction;
     const double rh_energyf = ((float)rh_rawenergy) * ((float)rh_fraction);
     if (!edm::isFinite(rh_energy)) {
       throw cms::Exception("PFClusterAlgo") << "rechit " << refhit->detId() << " has non-finite energy... "
                                             << "The input of the particle flow clustering seems to be corrupted.";
     }
     cl_energy += rh_energy;
     cl_energy_float += rh_energyf;
     // If time resolution is given, calculate weighted average
     if (_timeResolutionCalc) {
       const double res2 = 1. / _timeResolutionCalc->timeResolution2(rh_rawenergy);
       cl_time += rh_fraction * refhit->time() * res2;
       cl_timeweight += rh_fraction * res2;
     } else {  // assume resolution ~ 1/E**2
       const double rh_rawenergy2 = rh_rawenergy * rh_rawenergy;
       cl_timeweight += rh_rawenergy2 * rh_fraction;
       cl_time += rh_rawenergy2 * rh_fraction * refhit->time();
     }
     if (rh_energy > max_e) {
       max_e = rh_energy;
       max_e_layer = rhf.recHitRef()->layer();
       refmax = refhit;
     }
   }
   cluster.setEnergy(cl_energy);
   cluster.setTime(cl_time / cl_timeweight);
   if (_timeResolutionCalc) {
     cluster.setTimeError(std::sqrt(1.0f / float(cl_timeweight)));
   }
   cluster.setLayer(max_e_layer);
   const CaloSubdetectorGeometry* ecal_geom = nullptr;
   // get seed geometry information
   switch (max_e_layer) {
     case PFLayer::ECAL_BARREL:
       ecal_geom = _ebGeom;
       clusterT0 = _param_T0_EB;
       break;
     case PFLayer::ECAL_ENDCAP:
       ecal_geom = _eeGeom;
       clusterT0 = _param_T0_EE;
       break;
     default:
       throw cms::Exception("InvalidLayer") << "ECAL Position Calc only accepts ECAL_BARREL or ECAL_ENDCAP";
   }
 
   auto center_cell = ecal_geom->getGeometry(refmax->detId());
   const double ctreta = center_cell->etaPos();
   const double actreta = std::abs(ctreta);
   // need to change T0 if in ES
   if (actreta > preshowerStartEta && actreta < preshowerEndEta) {
     if (ctreta > 0 && _esPlus)
       clusterT0 = _param_T0_ES;
     if (ctreta < 0 && _esMinus)
       clusterT0 = _param_T0_ES;
   }
   // floats to reproduce exactly the EGM code
   const float maxDepth = _param_X0 * (clusterT0 + vdt::fast_log(cl_energy_float));
   const float maxToFront = center_cell->getPosition().mag();
   // calculate the position
   const double logETot_inv = -vdt::fast_log(cl_energy_float);
   double position_norm = 0.0;
   double x(0.0), y(0.0), z(0.0);
   for (const reco::PFRecHitFraction& rhf : cluster.recHitFractions()) {
     double weight = 0.0;
     const reco::PFRecHitRef& refhit = rhf.recHitRef();
     const double rh_energy = ((float)refhit->energy()) * ((float)rhf.fraction());
     if (rh_energy > 0.0)
       weight = std::max(0.0, (_param_W0 + vdt::fast_log(rh_energy) + logETot_inv));
     auto cell = ecal_geom->getGeometry(refhit->detId());
     const float depth = maxDepth + maxToFront - cell->getPosition().mag();
     const GlobalPoint pos = static_cast<const TruncatedPyramid*>(cell.get())->getPosition(depth);
 
     x += weight * pos.x();
     y += weight * pos.y();
     z += weight * pos.z();
 
     position_norm += weight;
   }
 
   // FALL BACK to LINEAR WEIGHTS
   if (position_norm == 0.) {
     for (const reco::PFRecHitFraction& rhf : cluster.recHitFractions()) {
       double weight = 0.0;
       const reco::PFRecHitRef& refhit = rhf.recHitRef();
       const double rh_energy = ((float)refhit->energy()) * ((float)rhf.fraction());
       if (rh_energy > 0.0)
         weight = rh_energy / cluster.energy();
 
       auto cell = ecal_geom->getGeometry(refhit->detId());
       const float depth = maxDepth + maxToFront - cell->getPosition().mag();
       const GlobalPoint pos = cell->getPosition(depth);
 
       x += weight * pos.x();
       y += weight * pos.y();
       z += weight * pos.z();
 
       position_norm += weight;
     }
   }
 
   if (position_norm < _minAllowedNorm) {
     edm::LogError("WeirdClusterNormalization") << "PFCluster too far from seeding cell: set position to (0,0,0).";
     cluster.setPosition(math::XYZPoint(0, 0, 0));
   } else {
     const double norm_inverse = 1.0 / position_norm;
     x *= norm_inverse;
     y *= norm_inverse;
     z *= norm_inverse;
 
     cluster.setPosition(math::XYZPoint(x, y, z));
     cluster.calculatePositionREP();
   }
 }

◆ calculateAndSetPositions()

void ECAL2DPositionCalcWithDepthCorr::calculateAndSetPositions	(	reco::PFClusterCollection &	clusters,
		const HcalPFCuts *
	)

overridevirtual

Implements PFCPositionCalculatorBase.

Definition at line 99 of file ECAL2DPositionCalcWithDepthCorr.cc.

References calculateAndSetPositionActual(), and bsc_activity_cfg::clusters.

                                                                                                                    {
   for (reco::PFCluster& cluster : clusters) {
     calculateAndSetPositionActual(cluster);
   }
 }

◆ operator=()

ECAL2DPositionCalcWithDepthCorr& ECAL2DPositionCalcWithDepthCorr::operator= ( const ECAL2DPositionCalcWithDepthCorr & )

delete

◆ update()

void ECAL2DPositionCalcWithDepthCorr::update ( const edm::EventSetup & es )

overridevirtual

Reimplemented from PFCPositionCalculatorBase.

Definition at line 80 of file ECAL2DPositionCalcWithDepthCorr.cc.

References _ebGeom, _eeGeom, _esGeom, _esMinus, _esPlus, _geomToken, DetId::Ecal, EcalBarrel, EcalEndcap, EcalPreshower, CaloSubdetectorGeometry::getGeometry(), edm::EventSetup::getHandle(), CaloGeometry::getSubdetectorGeometry(), CaloSubdetectorGeometry::getValidDetIds(), and z.

Referenced by progressbar.ProgressBar::__next__(), MatrixUtil.Matrix::__setitem__(), MatrixUtil.Steps::__setitem__(), progressbar.ProgressBar::finish(), and MatrixUtil.Steps::overwrite().

                                                                     {
   edm::ESHandle<CaloGeometry> geohandle = es.getHandle(_geomToken);
   _ebGeom = geohandle->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
   _eeGeom = geohandle->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
   _esGeom = geohandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
   if (_esGeom) {
     // ripped from RecoEcal/EgammaCoreTools
     for (uint32_t ic = 0; ic < _esGeom->getValidDetIds().size() && (!_esPlus || !_esMinus); ++ic) {
       const double z = _esGeom->getGeometry(_esGeom->getValidDetIds()[ic])->getPosition().z();
       _esPlus = _esPlus || (0 < z);
       _esMinus = _esMinus || (0 > z);
     }
   }
 }