#include <Basic2DGenericPFlowPositionCalc.h>

Inheritance diagram for Basic2DGenericPFlowPositionCalc:

Public Member Functions
	Basic2DGenericPFlowPositionCalc (const edm::ParameterSet &conf)

	Basic2DGenericPFlowPositionCalc (const Basic2DGenericPFlowPositionCalc &)=delete

void	calculateAndSetPosition (reco::PFCluster &) override

void	calculateAndSetPositions (reco::PFClusterCollection &) override

Basic2DGenericPFlowPositionCalc &	operator= (const Basic2DGenericPFlowPositionCalc &)=delete

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

PosCalc &	operator= (const PosCalc &)=delete

	PFCPositionCalculatorBase (const edm::ParameterSet &conf)

	PFCPositionCalculatorBase (const PosCalc &)=delete

virtual void	update (const edm::EventSetup &)

virtual	~PFCPositionCalculatorBase ()=default

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

Private Attributes
std::tuple< std::vector< int >, std::vector< int >, std::vector< float > >	_logWeightDenom

const float	_minAllowedNorm

const int	_posCalcNCrystals

std::unique_ptr< CaloRecHitResolutionProvider >	_timeResolutionCalcBarrel

std::unique_ptr< CaloRecHitResolutionProvider >	_timeResolutionCalcEndcap

Additional Inherited Members
Protected Attributes inherited from PFCPositionCalculatorBase
const float	_minFractionInCalc

Detailed Description

Definition at line 13 of file Basic2DGenericPFlowPositionCalc.h.

Constructor & Destructor Documentation

Basic2DGenericPFlowPositionCalc::Basic2DGenericPFlowPositionCalc ( const edm::ParameterSet & conf )

inline

Definition at line 15 of file Basic2DGenericPFlowPositionCalc.h.

References _logWeightDenom, _posCalcNCrystals, _timeResolutionCalcBarrel, _timeResolutionCalcEndcap, calculateAndSetPosition(), calculateAndSetPositions(), customizeHLTforCMSSW::depths, HLT_2018_cff::detectorEnum, Exception, edm::ParameterSet::exists(), edm::ParameterSet::getParameter(), edm::ParameterSet::getParameterSet(), edm::ParameterSet::getParameterSetVector(), mps_fire::i, operator=(), muonDTDigis_cfi::pset, and AlCaHLTBitMon_QueryRunRegistry::string.

       : PFCPositionCalculatorBase(conf),
         _posCalcNCrystals(conf.getParameter<int>("posCalcNCrystals")),
         _minAllowedNorm(conf.getParameter<double>("minAllowedNormalization")) {
     std::vector<int> detectorEnum;
     std::vector<int> depths;
     std::vector<double> logWeightDenom;
     std::vector<float> logWeightDenomInv;
 
     if (conf.exists("logWeightDenominatorByDetector")) {
       const std::vector<edm::ParameterSet>& logWeightDenominatorByDetectorPSet =
           conf.getParameterSetVector("logWeightDenominatorByDetector");
 
       for (const auto& pset : logWeightDenominatorByDetectorPSet) {
         if (!pset.exists("detector")) {
           throw cms::Exception("logWeightDenominatorByDetectorPSet") << "logWeightDenominator : detector not specified";
         }
 
         const std::string& det = pset.getParameter<std::string>("detector");
 
         if (det == std::string("HCAL_BARREL1") || det == std::string("HCAL_ENDCAP")) {
           std::vector<int> depthsT = pset.getParameter<std::vector<int> >("depths");
           std::vector<double> logWeightDenomT = pset.getParameter<std::vector<double> >("logWeightDenominator");
           if (logWeightDenomT.size() != depthsT.size()) {
             throw cms::Exception("logWeightDenominator") << "logWeightDenominator mismatch with the numbers of depths";
           }
           for (unsigned int i = 0; i < depthsT.size(); ++i) {
             if (det == std::string("HCAL_BARREL1"))
               detectorEnum.push_back(1);
             if (det == std::string("HCAL_ENDCAP"))
               detectorEnum.push_back(2);
             depths.push_back(depthsT[i]);
             logWeightDenom.push_back(logWeightDenomT[i]);
           }
         }
       }
     } else {
       detectorEnum.push_back(0);
       depths.push_back(0);
       logWeightDenom.push_back(conf.getParameter<double>("logWeightDenominator"));
     }
 
     for (unsigned int i = 0; i < depths.size(); ++i) {
       logWeightDenomInv.push_back(1. / logWeightDenom[i]);
     }
 
     //    _logWeightDenom = std::make_pair(depths,logWeightDenomInv);
     _logWeightDenom = std::make_tuple(detectorEnum, depths, logWeightDenomInv);
 
     _timeResolutionCalcBarrel.reset(nullptr);
     if (conf.exists("timeResolutionCalcBarrel")) {
       const edm::ParameterSet& timeResConf = conf.getParameterSet("timeResolutionCalcBarrel");
       _timeResolutionCalcBarrel.reset(new CaloRecHitResolutionProvider(timeResConf));
     }
     _timeResolutionCalcEndcap.reset(nullptr);
     if (conf.exists("timeResolutionCalcEndcap")) {
       const edm::ParameterSet& timeResConf = conf.getParameterSet("timeResolutionCalcEndcap");
       _timeResolutionCalcEndcap.reset(new CaloRecHitResolutionProvider(timeResConf));
     }
 
     switch (_posCalcNCrystals) {
       case 5:
       case 9:
       case -1:
         break;
       default:
         edm::LogError("Basic2DGenericPFlowPositionCalc") << "posCalcNCrystals not valid";
         assert(0);  // bug
     }
   }

Basic2DGenericPFlowPositionCalc::Basic2DGenericPFlowPositionCalc ( const Basic2DGenericPFlowPositionCalc & )

delete

Member Function Documentation

void Basic2DGenericPFlowPositionCalc::calculateAndSetPosition ( reco::PFCluster & cluster )

overridevirtual

Implements PFCPositionCalculatorBase.

Definition at line 20 of file Basic2DGenericPFlowPositionCalc.cc.

References calculateAndSetPositionActual().

Referenced by Basic2DGenericPFlowPositionCalc().

                                                                                     {
   calculateAndSetPositionActual(cluster);
 }

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

private

Definition at line 30 of file Basic2DGenericPFlowPositionCalc.cc.

References _logWeightDenom, _minAllowedNorm, PFCPositionCalculatorBase::_minFractionInCalc, _posCalcNCrystals, _timeResolutionCalcBarrel, _timeResolutionCalcEndcap, a, b, electrons_cff::bool, reco::PFCluster::calculatePositionREP(), bookConverter::compute(), declareDynArray, LEDCalibrationChannels::depth, reco::PFRecHit::depth(), HLT_2018_cff::detectorEnum, reco::PFRecHit::detId(), HCALHighEnergyHPDFilter_cfi::energy, Exception, f, myMath::fast_logf(), dqmMemoryStats::float, HLT_2018_cff::fraction, PFLayer::HCAL_BARREL1, PFLayer::HCAL_ENDCAP, photonIsolationHIProducer_cfi::hf, hfClusterShapes_cfi::hits, mps_fire::i, createfilelist::int, PixelPluginsPhase0_cfi::isBarrel, edm::isNotFinite(), dqmiolumiharvest::j, dqmdumpme::k, reco::PFRecHit::layer(), SiStripPI::max, ecalTB2006H4_GenSimDigiReco_cfg::mySeed, nhits, PFLayer::NONE, AlCaHLTBitMon_ParallelJobs::p, reco::PFRecHit::position(), reco::PFCluster::recHitFractions(), reco::CaloCluster::seed(), reco::PFCluster::setDepth(), reco::CaloCluster::setEnergy(), reco::PFCluster::setLayer(), reco::CaloCluster::setPosition(), reco::PFCluster::setTime(), reco::PFCluster::setTimeError(), mathSSE::sqrt(), MessageLogger_cff::threshold, reco::PFRecHit::time(), mitigatedMETSequence_cff::U, unInitDynArray, UNLIKELY, x, y, and z.

Referenced by calculateAndSetPosition(), and calculateAndSetPositions().

                                                                                                 {
   if (!cluster.seed()) {
     throw cms::Exception("ClusterWithNoSeed") << " Found a cluster with no seed: " << cluster;
   }
   double cl_energy = 0;
   double cl_time = 0;
   double cl_timeweight = 0.0;
   double max_e = 0.0;
   PFLayer::Layer max_e_layer = PFLayer::NONE;
   // find the seed and max layer and also calculate time
   //Michalis : Even if we dont use timing in clustering here we should fill
   //the time information for the cluster. This should use the timing resolution(1/E)
   //so the weight should be fraction*E^2
   //calculate a simplistic depth now. The log weighted will be done
   //in different stage
 
   auto const recHitCollection =
       &(*cluster.recHitFractions()[0].recHitRef()) - cluster.recHitFractions()[0].recHitRef().key();
   auto nhits = cluster.recHitFractions().size();
   struct LHit {
     reco::PFRecHit const* hit;
     float energy;
     float fraction;
   };
   declareDynArray(LHit, nhits, hits);
   for (auto i = 0U; i < nhits; ++i) {
     auto const& hf = cluster.recHitFractions()[i];
     auto k = hf.recHitRef().key();
     auto p = recHitCollection + k;
     hits[i] = {p, (*p).energy(), float(hf.fraction())};
   }
 
   bool resGiven = bool(_timeResolutionCalcBarrel) & bool(_timeResolutionCalcEndcap);
   LHit mySeed = {};
   for (auto const& rhf : hits) {
     const reco::PFRecHit& refhit = *rhf.hit;
     if (refhit.detId() == cluster.seed())
       mySeed = rhf;
     const auto rh_fraction = rhf.fraction;
     const auto rh_rawenergy = rhf.energy;
     const auto rh_energy = rh_rawenergy * rh_fraction;
 #ifdef PF_DEBUG
     if
       UNLIKELY(edm::isNotFinite(rh_energy)) {
         throw cms::Exception("PFClusterAlgo") << "rechit " << refhit.detId() << " has a NaN energy... "
                                               << "The input of the particle flow clustering seems to be corrupted.";
       }
 #endif
     cl_energy += rh_energy;
     // If time resolution is given, calculated weighted average
     if (resGiven) {
       double res2 = 1.e-4;
       int cell_layer = (int)refhit.layer();
       res2 = isBarrel(cell_layer) ? 1. / _timeResolutionCalcBarrel->timeResolution2(rh_rawenergy)
                                   : 1. / _timeResolutionCalcEndcap->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 = refhit.layer();
     }
   }
 
   cluster.setEnergy(cl_energy);
   cluster.setTime(cl_time / cl_timeweight);
   if (resGiven) {
     cluster.setTimeError(std::sqrt(1.0f / float(cl_timeweight)));
   }
   cluster.setLayer(max_e_layer);
 
   // calculate the position
   double depth = 0.0;
   double position_norm = 0.0;
   double x(0.0), y(0.0), z(0.0);
   if (nullptr != mySeed.hit) {
     auto seedNeighbours = mySeed.hit->neighbours();
     switch (_posCalcNCrystals) {
       case 5:
         seedNeighbours = mySeed.hit->neighbours4();
         break;
       case 9:
         seedNeighbours = mySeed.hit->neighbours8();
         break;
       default:
         break;
     }
 
     auto compute = [&](LHit const& rhf) {
       const reco::PFRecHit& refhit = *rhf.hit;
 
       int cell_layer = (int)refhit.layer();
       float threshold = 0;
 
       for (unsigned int j = 0; j < (std::get<2>(_logWeightDenom)).size(); ++j) {
         // barrel is detecor type1
         int detectorEnum = std::get<0>(_logWeightDenom)[j];
         int depth = std::get<1>(_logWeightDenom)[j];
 
         if ((cell_layer == PFLayer::HCAL_BARREL1 && detectorEnum == 1 && refhit.depth() == depth) ||
             (cell_layer == PFLayer::HCAL_ENDCAP && detectorEnum == 2 && refhit.depth() == depth) || detectorEnum == 0)
           threshold = std::get<2>(_logWeightDenom)[j];
       }
 
       const auto rh_energy = rhf.energy * rhf.fraction;
       const auto norm =
           (rhf.fraction < _minFractionInCalc ? 0.0f : std::max(0.0f, vdt::fast_logf(rh_energy * threshold)));
       const auto rhpos_xyz = refhit.position() * norm;
       x += rhpos_xyz.x();
       y += rhpos_xyz.y();
       z += rhpos_xyz.z();
       depth += refhit.depth() * norm;
       position_norm += norm;
     };
 
     if (_posCalcNCrystals != -1)  // sorted to make neighbour search faster (maybe)
       std::sort(hits.begin(), hits.end(), [](LHit const& a, LHit const& b) { return a.hit < b.hit; });
 
     if (_posCalcNCrystals == -1)
       for (auto const& rhf : hits)
         compute(rhf);
     else {  // only seed and its neighbours
       compute(mySeed);
       // search seedNeighbours to find energy fraction in cluster (sic)
       unInitDynArray(reco::PFRecHit const*, seedNeighbours.size(), nei);
       for (auto k : seedNeighbours) {
         nei.push_back(recHitCollection + k);
       }
       std::sort(nei.begin(), nei.end());
       struct LHitLess {
         auto operator()(LHit const& a, reco::PFRecHit const* b) const { return a.hit < b; }
         auto operator()(reco::PFRecHit const* b, LHit const& a) const { return b < a.hit; }
       };
       std::set_intersection(
           hits.begin(), hits.end(), nei.begin(), nei.end(), boost::make_function_output_iterator(compute), LHitLess());
     }
   } else {
     throw cms::Exception("Basic2DGenerticPFlowPositionCalc")
         << "Cluster seed hit is null, something is wrong with PFlow RecHit!";
   }
 
   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));
     cluster.calculatePositionREP();
   } else {
     const double norm_inverse = 1.0 / position_norm;
     x *= norm_inverse;
     y *= norm_inverse;
     z *= norm_inverse;
     depth *= norm_inverse;
     cluster.setPosition(math::XYZPoint(x, y, z));
     cluster.setDepth(depth);
     cluster.calculatePositionREP();
   }
 }

void Basic2DGenericPFlowPositionCalc::calculateAndSetPositions ( reco::PFClusterCollection & clusters )

overridevirtual

Implements PFCPositionCalculatorBase.

Definition at line 24 of file Basic2DGenericPFlowPositionCalc.cc.

References calculateAndSetPositionActual().

Referenced by Basic2DGenericPFlowPositionCalc().

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