#include <Basic2DGenericPFlowPositionCalc.h>

Inheritance diagram for Basic2DGenericPFlowPositionCalc:

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

	Basic2DGenericPFlowPositionCalc (const Basic2DGenericPFlowPositionCalc &)=delete

void	calculateAndSetPosition (reco::PFCluster &)

void	calculateAndSetPositions (reco::PFClusterCollection &)

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 &)

	~PFCPositionCalculatorBase ()

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

Private Attributes
const 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 12 of file Basic2DGenericPFlowPositionCalc.h.

Constructor & Destructor Documentation

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

inline

Definition at line 14 of file Basic2DGenericPFlowPositionCalc.h.

References _posCalcNCrystals, _timeResolutionCalcBarrel, _timeResolutionCalcEndcap, assert(), edm::ParameterSet::exists(), edm::ParameterSet::getParameterSet(), and NULL.

                                                                :
     PFCPositionCalculatorBase(conf),    
     _posCalcNCrystals(conf.getParameter<int>("posCalcNCrystals")),
     _logWeightDenom(1./conf.getParameter<double>("logWeightDenominator")),
     _minAllowedNorm(conf.getParameter<double>("minAllowedNormalization"))
 
   {  
     _timeResolutionCalcBarrel.reset(NULL);
     if( conf.exists("timeResolutionCalcBarrel") ) {
       const edm::ParameterSet& timeResConf = 
         conf.getParameterSet("timeResolutionCalcBarrel");
         _timeResolutionCalcBarrel.reset(new CaloRecHitResolutionProvider(timeResConf));
     }
     _timeResolutionCalcEndcap.reset(NULL);
     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 )

virtual

Implements PFCPositionCalculatorBase.

Definition at line 24 of file Basic2DGenericPFlowPositionCalc.cc.

References calculateAndSetPositionActual().

                                                 {
   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, reco::PFCluster::calculatePositionREP(), bookConverter::compute(), declareDynArray, reco::PFRecHit::depth(), HLT_25ns10e33_v2_cff::depth, reco::PFRecHit::detId(), relval_parameters_module::energy, reco::PFRecHit::energy(), Exception, f, myMath::fast_logf(), HLT_25ns10e33_v2_cff::fraction, i, GeomDetEnumerators::isBarrel(), edm::isNotFinite(), relval_2017::k, reco::PFRecHit::layer(), bookConverter::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::PFRecHit::time(), 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);
   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);
   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;  
     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*_logWeightDenom)) );
     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()
      );
   }
 
 
   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 )

virtual

Implements PFCPositionCalculatorBase.

Definition at line 29 of file Basic2DGenericPFlowPositionCalc.cc.

References calculateAndSetPositionActual().

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