#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(), particleFlowClusterHBHE_cfi::depths, particleFlowRecHitHBHE_cfi::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 24 of file Basic2DGenericPFlowPositionCalc.cc.

References calculateAndSetPositionActual(), and calculateAndSetPositions().

Referenced by Basic2DGenericPFlowPositionCalc().

                                                 {
   calculateAndSetPositionActual(cluster);
 }

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

private

Definition at line 36 of file Basic2DGenericPFlowPositionCalc.cc.

References _logWeightDenom, _minAllowedNorm, PFCPositionCalculatorBase::_minFractionInCalc, _posCalcNCrystals, _timeResolutionCalcBarrel, _timeResolutionCalcEndcap, a, b, electrons_cff::bool, reco::PFCluster::calculatePositionREP(), bookConverter::compute(), declareDynArray, egammaForCoreTracking_cff::depth, reco::PFRecHit::depth(), particleFlowRecHitHBHE_cfi::detectorEnum, reco::PFRecHit::detId(), Exception, f, myMath::fast_logf(), objects.autophobj::float, dedxEstimators_cff::fraction, PFLayer::HCAL_BARREL1, PFLayer::HCAL_ENDCAP, photonIsolationHIProducer_cfi::hf, hfClusterShapes_cfi::hits, mps_fire::i, createfilelist::int, GeomDetEnumerators::isBarrel(), edm::isNotFinite(), gen::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(), jetUpdater_cfi::sort, mathSSE::sqrt(), electronIdCutBased_cfi::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={nullptr};
   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 29 of file Basic2DGenericPFlowPositionCalc.cc.

References calculateAndSetPositionActual().

Referenced by Basic2DGenericPFlowPositionCalc(), and calculateAndSetPosition().

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