An explicit threshold-based clustering algorithm. More...

#include <MTDThresholdClusterizer.h>

Inheritance diagram for MTDThresholdClusterizer:

Public Member Functions
void	clusterize (const FTLRecHitCollection &input, const MTDGeometry geom, const MTDTopology topo, FTLClusterCollection &output) override
	Cluster hits. This method operates on a matrix of hits and finds the largest contiguous cluster around each seed hit. Input and output data stored in DetSet. More...

	MTDThresholdClusterizer (edm::ParameterSet const &conf)
	Constructor: More...

	~MTDThresholdClusterizer () override

Public Member Functions inherited from MTDClusterizerBase
virtual	~MTDClusterizerBase ()

Static Public Member Functions
static void	fillPSetDescription (edm::ParameterSetDescription &desc)

Private Member Functions
void	clear_buffer (RecHitIterator itr)
	Clear the internal buffer array. More...

void	copy_to_buffer (RecHitIterator itr, const MTDGeometry geom, const MTDTopology topo)
	Copy FTLRecHit into the buffer, identify seeds. More...

FTLCluster	make_cluster (const FTLCluster::FTLHitPos &hit)
	The actual clustering algorithm: group the neighboring hits around the seed. More...

bool	setup (const MTDGeometry geometry, const MTDTopology topo, const DetId &id)

Private Attributes
bool	bufferAlreadySet

MTDArrayBuffer	theBuffer
	Data storage. More...

std::vector< FTLCluster >	theClusters

const float	theClusterThreshold

DetId	theCurrentId

const float	theHitThreshold
	Clustering-related quantities: More...

int	theNumOfCols

int	theNumOfRows
	Geometry-related information. More...

const float	thePositionThreshold

std::vector< FTLCluster::FTLHitPos >	theSeeds

const float	theSeedThreshold

const float	theTimeThreshold

Additional Inherited Members
Public Types inherited from MTDClusterizerBase
typedef FTLClusterCollection::const_iterator	ClusterIterator

typedef FTLRecHitCollection::const_iterator	RecHitIterator

Detailed Description

An explicit threshold-based clustering algorithm.

A threshold-based clustering algorithm which clusters FTLRecHits into FTLClusters for each DetUnit. The algorithm is straightforward and purely topological: the clustering process starts with seed hits and continues by adding adjacent hits above the hit threshold. Once the cluster is made, it has to be above the cluster threshold as well.

The clusterization is performed on a matrix with size equal to the size of the MTD detector, each cell containing the cahrge and time of the corresponding hit The matrix is reset after each clusterization.

The search starts from seed hits, i.e. hits with sufficiently large amplitudes

FTLCluster contains a barycenter, but it should be noted that that information is largely useless. One must use a PositionEstimator class to compute the RecHit position and its error for every given cluster.

Sets the MTDArrayBuffer dimensions and pixel thresholds. Makes clusters and stores them in theCache if the option useCache has been set.

Definition at line 51 of file MTDThresholdClusterizer.h.

Constructor & Destructor Documentation

◆ MTDThresholdClusterizer()

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

Constructor:

Definition at line 26 of file MTDThresholdClusterizer.cc.

     :  // Get energy thresholds
       theHitThreshold(conf.getParameter<double>("HitThreshold")),
       theSeedThreshold(conf.getParameter<double>("SeedThreshold")),
       theClusterThreshold(conf.getParameter<double>("ClusterThreshold")),
       theTimeThreshold(conf.getParameter<double>("TimeThreshold")),
       thePositionThreshold(conf.getParameter<double>("PositionThreshold")),
       theNumOfRows(0),
       theNumOfCols(0),
       theCurrentId(0),
       theBuffer(theNumOfRows, theNumOfCols),
       bufferAlreadySet(false) {}

◆ ~MTDThresholdClusterizer()

MTDThresholdClusterizer::~MTDThresholdClusterizer ( )

override

Definition at line 39 of file MTDThresholdClusterizer.cc.

39 {}

Member Function Documentation

◆ clear_buffer()

void MTDThresholdClusterizer::clear_buffer ( RecHitIterator itr )

private

Clear the internal buffer array.

MTDs which are not part of recognized clusters are NOT ERASED during the cluster finding. Erase them now.

Definition at line 189 of file MTDThresholdClusterizer.cc.

References MTDArrayBuffer::clear(), and theBuffer.

Referenced by clusterize().

189 { theBuffer.clear(itr->row(), itr->column()); }

MTDThresholdClusterizer::theBuffer

MTDArrayBuffer theBuffer

Data storage.

Definition: MTDThresholdClusterizer.h:82

MTDArrayBuffer::clear

void clear(uint row, uint col)

Definition: MTDArrayBuffer.h:52

◆ clusterize()

void MTDThresholdClusterizer::clusterize	(	const FTLRecHitCollection &	input,
		const MTDGeometry *	geom,
		const MTDTopology *	topo,
		FTLClusterCollection &	output
	)

overridevirtual

Cluster hits. This method operates on a matrix of hits and finds the largest contiguous cluster around each seed hit. Input and output data stored in DetSet.

Implements MTDClusterizerBase.

Definition at line 90 of file MTDThresholdClusterizer.cc.

References cms::cuda::assert(), MTDDetId::BTL, clear_buffer(), copy_to_buffer(), MTDTopologyMode::crysLayoutFromTopoMode(), TauDecayModes::dec, mps_fire::end, MTDArrayBuffer::energy(), FTLCluster::energy(), MTDDetId::ETL, Exception, MTDDetId::FastTime, relativeConstraints::geom, FTLCluster::getClusterErrorX(), FTLCluster::getClusterPosX(), MTDTopology::getMTDTopologyMode(), mps_fire::i, input, LogDebug, make_cluster(), MTDDetId::mtdSubDetector(), edmNew::DetSetVector< T >::FastFiller::push_back(), FastTimerService_cff::range, DetId::rawId(), setup(), FTLCluster::size(), MTDDetId::subDetector(), theBuffer, theClusterThreshold, theSeeds, theSeedThreshold, FTLCluster::time(), FTLCluster::timeError(), FTLCluster::x(), and FTLCluster::y().

                                                                        {
   FTLRecHitCollection::const_iterator begin = input.begin();
   FTLRecHitCollection::const_iterator end = input.end();
 
   // Do not bother for empty detectors
   if (begin == end) {
     edm::LogInfo("MTDThresholdClusterizer") << "No hits to clusterize";
     return;
   }
 
   LogDebug("MTDThresholdClusterizer") << "Input collection " << input.size();
   assert(output.empty());
 
   std::set<unsigned> geoIds;
   std::multimap<unsigned, unsigned> geoIdToIdx;
 
   unsigned index = 0;
   for (const auto& hit : input) {
     MTDDetId mtdId = MTDDetId(hit.detid());
     if (mtdId.subDetector() != MTDDetId::FastTime) {
       throw cms::Exception("MTDThresholdClusterizer")
           << "MTDDetId: " << std::hex << mtdId.rawId() << " is invalid!" << std::dec << std::endl;
     }
 
     if (mtdId.mtdSubDetector() == MTDDetId::BTL) {
       BTLDetId hitId(hit.detid());
       //for BTL topology gives different layout id
       DetId geoId = hitId.geographicalId(MTDTopologyMode::crysLayoutFromTopoMode(topo->getMTDTopologyMode()));
       geoIdToIdx.emplace(geoId, index);
       geoIds.emplace(geoId);
       ++index;
     } else if (mtdId.mtdSubDetector() == MTDDetId::ETL) {
       ETLDetId hitId(hit.detid());
       DetId geoId = hitId.geographicalId();
       geoIdToIdx.emplace(geoId, index);
       geoIds.emplace(geoId);
       ++index;
     } else {
       throw cms::Exception("MTDThresholdClusterizer")
           << "MTDDetId: " << std::hex << mtdId.rawId() << " is invalid!" << std::dec << std::endl;
     }
   }
 
   //cluster hits within geoIds (modules)
   for (unsigned id : geoIds) {
     //  Set up the clusterization on this DetId.
     if (!setup(geom, topo, DetId(id)))
       return;
 
     auto range = geoIdToIdx.equal_range(id);
     LogDebug("MTDThresholdClusterizer") << "Matching Ids for " << std::hex << id << std::dec << " ["
                                         << range.first->second << "," << range.second->second << "]";
 
     //  Copy MTDRecHits to the buffer array; select the seed hits
     //  on the way, and store them in theSeeds.
     for (auto itr = range.first; itr != range.second; ++itr) {
       const unsigned hitidx = itr->second;
       copy_to_buffer(begin + hitidx, geom, topo);
     }
 
     FTLClusterCollection::FastFiller clustersOnDet(output, id);
 
     for (unsigned int i = 0; i < theSeeds.size(); i++) {
       if (theBuffer.energy(theSeeds[i]) > theSeedThreshold) {  // Is this seed still valid?
         //  Make a cluster around this seed
         const FTLCluster& cluster = make_cluster(theSeeds[i]);
 
         //  Check if the cluster is above threshold
         if (cluster.energy() > theClusterThreshold) {
           LogDebug("MTDThresholdClusterizer")
               << "putting in this cluster " << i << " #hits:" << cluster.size() << " E:" << cluster.energy()
               << " T +/- DT:" << cluster.time() << " +/- " << cluster.timeError() << " X:" << cluster.x()
               << " Y:" << cluster.y() << " xpos +/- err " << cluster.getClusterPosX() << " +/- "
               << cluster.getClusterErrorX();
           clustersOnDet.push_back(cluster);
         }
       }
     }
 
     // Erase the seeds.
     theSeeds.clear();
     //  Need to clean unused hits from the buffer array.
     for (auto itr = range.first; itr != range.second; ++itr) {
       const unsigned hitidx = itr->second;
       clear_buffer(begin + hitidx);
     }
   }
 }

◆ copy_to_buffer()

void MTDThresholdClusterizer::copy_to_buffer	(	RecHitIterator	itr,
		const MTDGeometry *	geom,
		const MTDTopology *	topo
	)

private

Copy FTLRecHit into the buffer, identify seeds.

Definition at line 194 of file MTDThresholdClusterizer.cc.

References GeomDetEnumerators::barrel, MTDDetId::BTL, cuy::col, MTDTopologyMode::crysLayoutFromTopoMode(), GeomDetEnumerators::endcap, hcalRecHitTable_cff::energy, MTDDetId::ETL, relativeConstraints::geom, MTDTopology::getMTDTopologyMode(), GeomDetEnumerators::invalidLoc, BTLRecHitsErrorEstimatorIM::localError(), LogDebug, MTDDetId::mtdSubDetector(), RectangularMTDTopology::pixelToModuleLocalPoint(), position, DetId::rawId(), MTDArrayBuffer::set(), ProxyMTDTopology::specificTopology(), mathSSE::sqrt(), theBuffer, theHitThreshold, theSeeds, theSeedThreshold, hcalRecHitTable_cff::time, PV3DBase< T, PVType, FrameType >::x(), LocalError::xx(), PV3DBase< T, PVType, FrameType >::y(), LocalError::yy(), and PV3DBase< T, PVType, FrameType >::z().

Referenced by clusterize().

                                                                                                                  {
   MTDDetId mtdId = MTDDetId(itr->detid());
   int row = itr->row();
   int col = itr->column();
   GeomDetEnumerators::Location subDet = GeomDetEnumerators::invalidLoc;
   float energy = itr->energy();
   float time = itr->time();
   float timeError = itr->timeError();
   float position = itr->position();
   float xpos = 0.f;
   // position is the longitudinal offset that should be added into local x for bars in phi geometry
   LocalPoint local_point(0, 0, 0);
   LocalError local_error(0, 0, 0);
   if (mtdId.mtdSubDetector() == MTDDetId::BTL) {
     subDet = GeomDetEnumerators::barrel;
     BTLDetId id = itr->id();
     DetId geoId = id.geographicalId(MTDTopologyMode::crysLayoutFromTopoMode(topo->getMTDTopologyMode()));
     const auto& det = geom->idToDet(geoId);
     const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(det->topology());
     const RectangularMTDTopology& topol = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
 
     LocalPoint lp_pixel(position, 0, 0);
     local_point = topol.pixelToModuleLocalPoint(lp_pixel, row, col);
     BTLRecHitsErrorEstimatorIM btlError(det, local_point);
     local_error = btlError.localError();
     xpos = local_point.x();
   } else if (mtdId.mtdSubDetector() == MTDDetId::ETL) {
     subDet = GeomDetEnumerators::endcap;
     ETLDetId id = itr->id();
     DetId geoId = id.geographicalId();
     const auto& det = geom->idToDet(geoId);
     const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(det->topology());
     const RectangularMTDTopology& topol = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
 
     LocalPoint lp_pixel(0, 0, 0);
     local_point = topol.pixelToModuleLocalPoint(lp_pixel, row, col);
   }
 
   LogDebug("MTDThresholdClusterizer") << "DetId " << mtdId.rawId() << " subd " << mtdId.mtdSubDetector() << " row/col "
                                       << row << " / " << col << " energy " << energy << " time " << time
                                       << " time error " << timeError << " local_point " << local_point.x() << " "
                                       << local_point.y() << " " << local_point.z() << " local error "
                                       << std::sqrt(local_error.xx()) << " " << std::sqrt(local_error.yy()) << " xpos "
                                       << xpos;
   if (energy > theHitThreshold) {
     theBuffer.set(row, col, subDet, energy, time, timeError, local_error, local_point, xpos);
     if (energy > theSeedThreshold)
       theSeeds.push_back(FTLCluster::FTLHitPos(row, col));
     //sort seeds?
   }
 }

◆ fillPSetDescription()

void MTDThresholdClusterizer::fillPSetDescription ( edm::ParameterSetDescription & desc )

static

Definition at line 42 of file MTDThresholdClusterizer.cc.

References submitPVResolutionJobs::desc.

Referenced by MTDClusterProducer::fillDescriptions().

                                                                                   {
   desc.add<double>("HitThreshold", 0.);
   desc.add<double>("SeedThreshold", 0.);
   desc.add<double>("ClusterThreshold", 0.);
   desc.add<double>("TimeThreshold", 10.);
   desc.add<double>("PositionThreshold", -1.0);
 }

◆ make_cluster()

FTLCluster MTDThresholdClusterizer::make_cluster ( const FTLCluster::FTLHitPos & hit )

private

The actual clustering algorithm: group the neighboring hits around the seed.

Definition at line 249 of file MTDThresholdClusterizer.cc.

Referenced by clusterize().

                                                                                {
   //First we acquire the seeds for the clusters
 
   GeomDetEnumerators::Location seed_subdet = theBuffer.subDet(hit.row(), hit.col());
   float seed_energy = theBuffer.energy(hit.row(), hit.col());
   float seed_time = theBuffer.time(hit.row(), hit.col());
   float seed_time_error = theBuffer.time_error(hit.row(), hit.col());
   auto const seedPoint = theBuffer.local_point(hit.row(), hit.col());
   double seed_error_xx = theBuffer.local_error(hit.row(), hit.col()).xx();
   double seed_error_yy = theBuffer.local_error(hit.row(), hit.col()).yy();
   double seed_xpos = theBuffer.xpos(hit.row(), hit.col());
   theBuffer.clear(hit);
 
   AccretionCluster acluster;
   acluster.add(hit, seed_energy, seed_time, seed_time_error);
 
   // for BTL position along crystals add auxiliary vectors
   std::array<float, AccretionCluster::MAXSIZE> pixel_x{{-99999.}};
   std::array<float, AccretionCluster::MAXSIZE> pixel_errx2{{-99999.}};
   if (seed_subdet == GeomDetEnumerators::barrel) {
     pixel_x[acluster.top()] = seed_xpos;
     pixel_errx2[acluster.top()] = seed_error_xx;
   }
 
   bool stopClus = false;
   //Here we search all hits adjacent to all hits in the cluster.
   while (!acluster.empty() && !stopClus) {
     //This is the standard algorithm to find and add a hit
     auto curInd = acluster.top();
     acluster.pop();
     for (auto c = std::max(0, int(acluster.y[curInd] - 1));
          c < std::min(int(acluster.y[curInd] + 2), int(theBuffer.columns())) && !stopClus;
          ++c) {
       for (auto r = std::max(0, int(acluster.x[curInd] - 1));
            r < std::min(int(acluster.x[curInd] + 2), int(theBuffer.rows())) && !stopClus;
            ++r) {
         LogDebug("MTDThresholdClusterizer")
             << "Clustering subdet " << seed_subdet << " around " << curInd << " X,Y = " << acluster.x[curInd] << ","
             << acluster.y[curInd] << " r,c = " << r << "," << c << " energy,time = " << theBuffer.energy(r, c) << " "
             << theBuffer.time(r, c);
         if (theBuffer.energy(r, c) > theHitThreshold) {
           if (std::abs(theBuffer.time(r, c) - seed_time) >
               theTimeThreshold *
                   sqrt(theBuffer.time_error(r, c) * theBuffer.time_error(r, c) + seed_time_error * seed_time_error))
             continue;
           if ((seed_subdet == GeomDetEnumerators::barrel) && (theBuffer.subDet(r, c) == GeomDetEnumerators::barrel) &&
               (thePositionThreshold > 0)) {
             double hit_error_xx = theBuffer.local_error(r, c).xx();
             double hit_error_yy = theBuffer.local_error(r, c).yy();
             if (((theBuffer.local_point(r, c) - seedPoint).mag2()) >
                 thePositionThreshold * thePositionThreshold *
                     (hit_error_xx + seed_error_xx + hit_error_yy + seed_error_yy))
               continue;
           }
           FTLCluster::FTLHitPos newhit(r, c);
           if (!acluster.add(newhit, theBuffer.energy(r, c), theBuffer.time(r, c), theBuffer.time_error(r, c))) {
             stopClus = true;
             break;
           }
           if (theBuffer.subDet(r, c) == GeomDetEnumerators::barrel) {
             pixel_x[acluster.top()] = theBuffer.xpos(r, c);
             pixel_errx2[acluster.top()] = theBuffer.local_error(r, c).xx();
           }
           theBuffer.clear(newhit);
         }
       }
     }
   }  // while accretion
 
   FTLCluster cluster(theCurrentId,
                      acluster.isize,
                      acluster.energy.data(),
                      acluster.time.data(),
                      acluster.timeError.data(),
                      acluster.x.data(),
                      acluster.y.data(),
                      acluster.xmin,
                      acluster.ymin);
 
   // For BTL compute the optimal position along crystal and uncertainty on it in absolute length units
 
   if (seed_subdet == GeomDetEnumerators::barrel) {
     float sumW(0.f), sumXW(0.f), sumXW2(0.f);
     for (unsigned int index = 0; index < acluster.top(); index++) {
       sumW += acluster.energy[index];
       sumXW += acluster.energy[index] * pixel_x[index];
       sumXW2 += acluster.energy[index] * acluster.energy[index] * pixel_errx2[index];
     }
     cluster.setClusterPosX(sumXW / sumW);
     cluster.setClusterErrorX(std::sqrt(sumXW2) / sumW);
   }
 
   return cluster;
 }

◆ setup()

bool MTDThresholdClusterizer::setup	(	const MTDGeometry *	geom,
		const MTDTopology *	topo,
		const DetId &	id
	)

private

Prepare the Clusterizer to work on a particular DetUnit. Re-init the size of the panel/plaquette (so update nrows and ncols),

Definition at line 54 of file MTDThresholdClusterizer.cc.

References bufferAlreadySet, MTDArrayBuffer::columns(), TauDecayModes::dec, Exception, relativeConstraints::geom, l1ctLayer2EG_cff::id, LogDebug, hgcalPlots::ncols, RectangularMTDTopology::ncolumns(), RectangularMTDTopology::nrows(), MTDArrayBuffer::rows(), MTDArrayBuffer::setSize(), ProxyMTDTopology::specificTopology(), theBuffer, theCurrentId, theNumOfCols, and theNumOfRows.

Referenced by clusterize().

                                                                                                      {
   theCurrentId = id;
   //using geopraphicalId here
   const auto& thedet = geom->idToDet(id);
   if (thedet == nullptr) {
     throw cms::Exception("MTDThresholdClusterizer")
         << "GeographicalID: " << std::hex << id.rawId() << " is invalid!" << std::dec << std::endl;
     return false;
   }
   const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(thedet->topology());
   const RectangularMTDTopology& topol = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
 
   // Get the new sizes.
   unsigned int nrows = topol.nrows();     // rows in x
   unsigned int ncols = topol.ncolumns();  // cols in y
 
   theNumOfRows = nrows;  // Set new sizes
   theNumOfCols = ncols;
 
   LogDebug("MTDThresholdClusterizer") << "Buffer size [" << theNumOfRows + 1 << "," << theNumOfCols + 1 << "]";
 
   if (nrows > theBuffer.rows() || ncols > theBuffer.columns()) {  // change only when a larger is needed
     // Resize the buffer
     theBuffer.setSize(nrows, ncols);
     bufferAlreadySet = true;
   }
 
   return true;
 }

Member Data Documentation

◆ bufferAlreadySet

bool MTDThresholdClusterizer::bufferAlreadySet

private

Definition at line 83 of file MTDThresholdClusterizer.h.

Referenced by setup().

◆ theBuffer

MTDArrayBuffer MTDThresholdClusterizer::theBuffer

private

Data storage.

Definition at line 82 of file MTDThresholdClusterizer.h.

Referenced by clear_buffer(), clusterize(), copy_to_buffer(), make_cluster(), and setup().

◆ theClusters

std::vector<FTLCluster> MTDThresholdClusterizer::theClusters

private

Definition at line 66 of file MTDThresholdClusterizer.h.

◆ theClusterThreshold

const float MTDThresholdClusterizer::theClusterThreshold

private

Definition at line 71 of file MTDThresholdClusterizer.h.

Referenced by clusterize().

◆ theCurrentId

DetId MTDThresholdClusterizer::theCurrentId

private

Definition at line 79 of file MTDThresholdClusterizer.h.

Referenced by make_cluster(), and setup().

◆ theHitThreshold

const float MTDThresholdClusterizer::theHitThreshold

private

Clustering-related quantities:

Definition at line 69 of file MTDThresholdClusterizer.h.

Referenced by copy_to_buffer(), and make_cluster().

◆ theNumOfCols

int MTDThresholdClusterizer::theNumOfCols

private

Definition at line 77 of file MTDThresholdClusterizer.h.

Referenced by setup().

◆ theNumOfRows

int MTDThresholdClusterizer::theNumOfRows

private

Geometry-related information.

Definition at line 76 of file MTDThresholdClusterizer.h.

Referenced by setup().

◆ thePositionThreshold

const float MTDThresholdClusterizer::thePositionThreshold

private

Definition at line 73 of file MTDThresholdClusterizer.h.

Referenced by make_cluster().

◆ theSeeds

std::vector<FTLCluster::FTLHitPos> MTDThresholdClusterizer::theSeeds

private

Definition at line 65 of file MTDThresholdClusterizer.h.

Referenced by clusterize(), and copy_to_buffer().

◆ theSeedThreshold

const float MTDThresholdClusterizer::theSeedThreshold

private

Definition at line 70 of file MTDThresholdClusterizer.h.

Referenced by clusterize(), and copy_to_buffer().

◆ theTimeThreshold

const float MTDThresholdClusterizer::theTimeThreshold

private

Definition at line 72 of file MTDThresholdClusterizer.h.

Referenced by make_cluster().

Public Member Functions

Static Public Member Functions

Private Member Functions

Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ MTDThresholdClusterizer()

◆ ~MTDThresholdClusterizer()

Member Function Documentation

◆ clear_buffer()

◆ clusterize()

◆ copy_to_buffer()

◆ fillPSetDescription()

◆ make_cluster()

◆ setup()

Member Data Documentation

◆ bufferAlreadySet

◆ theBuffer

◆ theClusters

◆ theClusterThreshold

◆ theCurrentId

◆ theHitThreshold

◆ theNumOfCols

◆ theNumOfRows

◆ thePositionThreshold

◆ theSeeds

◆ theSeedThreshold

◆ theTimeThreshold