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

#include <PixelThresholdClusterizer.h>

Inheritance diagram for PixelThresholdClusterizer:

Public Member Functions
void	clusterizeDetUnit (const edm::DetSet< PixelDigi > &input, const PixelGeomDetUnit *pixDet, const std::vector< short > &badChannels, edmNew::DetSetVector< SiPixelCluster >::FastFiller &output)
	Cluster pixels. This method operates on a matrix of pixels and finds the largest contiguous cluster around each seed pixel. Input and output data stored in DetSet.
	PixelThresholdClusterizer (edm::ParameterSet const &conf)
	~PixelThresholdClusterizer ()
Private Member Functions
int	calibrate (int adc, int col, int row)
void	clear_buffer (DigiIterator begin, DigiIterator end)
	Clear the internal buffer array.
void	copy_to_buffer (DigiIterator begin, DigiIterator end)
	Copy adc counts from PixelDigis into the buffer, identify seeds.
SiPixelCluster	make_cluster (const SiPixelCluster::PixelPos &pix, edmNew::DetSetVector< SiPixelCluster >::FastFiller &output)
	The actual clustering algorithm: group the neighboring pixels around the seed.
bool	setup (const PixelGeomDetUnit *pixDet)
	Private helper methods:
Private Attributes
bool	bufferAlreadySet
edm::ParameterSet	conf_
bool	dead_flag
uint32_t	detid_
bool	doMissCalibrate
bool	doSplitClusters
SiPixelArrayBuffer	theBuffer
	Data storage.
std::vector< SiPixelCluster >	theClusters
float	theClusterThreshold
float	theClusterThresholdInNoiseUnits
int	theConversionFactor
int	theFirstStack_
int	theNumOfCols
int	theNumOfRows
	Geometry-related information.
int	theOffset
int	thePixelThreshold
float	thePixelThresholdInNoiseUnits
	Clustering-related quantities:
std::vector < SiPixelCluster::PixelPos >	theSeeds
int	theSeedThreshold
float	theSeedThresholdInNoiseUnits
int	theStackADC_

Detailed Description

An explicit threshold-based clustering algorithm.

A specific threshold-based pixel clustering algorithm.

A threshold-based clustering algorithm which clusters SiPixelDigis into SiPixelClusters for each DetUnit. The algorithm is straightforward and purely topological: the clustering process starts with seed pixels and continues by adding adjacent pixels above the pixel 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 pixel detector, each cell containing the ADC count of the corresponding pixel. The matrix is reset after each clusterization.

The search starts from seed pixels, i.e. pixels with sufficiently large amplitudes, found at the time of filling of the matrix and stored in a

At this point the noise and dead channels are ignored, but soon they won't be.

SiPixelCluster contains a barrycenter, 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.

Author:: Largely copied from NewPixelClusterizer in ORCA written by Danek Kotlinski (PSI). Ported to CMSSW by Petar Maksimovic (JHU). DetSetVector data container implemented by V.Chiochia (Uni Zurich)

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

General logic of PixelThresholdClusterizer:

The clusterization is performed on a matrix with size equal to the size of the pixel detector, each cell containing the ADC count of the corresponding pixel. The matrix is reset after each clusterization.

The search starts from seed pixels, i.e. pixels with sufficiently large amplitudes, found at the time of filling of the matrix and stored in a SiPixelArrayBuffer.

Translate the pixel charge to electrons, we are suppose to do the calibrations ADC->electrons here. Modify the thresholds to be in electrons, convert adc to electrons. d.k. 20/3/06 Get rid of the noiseVector. d.k. 28/3/06

Definition at line 56 of file PixelThresholdClusterizer.h.

Constructor & Destructor Documentation

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

Constructor: Initilize the buffer to hold pixels from a detector module. This is a vector of 44k ints, stays valid all the time.

Definition at line 44 of file PixelThresholdClusterizer.cc.

References edm::ParameterSet::getParameter().

                                :
    conf_(conf), bufferAlreadySet(false), theNumOfRows(0), theNumOfCols(0), detid_(0) 
{
  // Get thresholds in electrons
  thePixelThreshold   = 
    conf_.getParameter<int>("ChannelThreshold");
  theSeedThreshold    = 
    conf_.getParameter<int>("SeedThreshold");
  theClusterThreshold = 
    conf_.getParameter<double>("ClusterThreshold");
  theConversionFactor = 
    conf_.getParameter<int>("VCaltoElectronGain");
  theOffset = 
    conf_.getParameter<int>("VCaltoElectronOffset");
  if ( conf_.exists("AdcFullScaleStack") ) theStackADC_=conf_.getParameter<int>("AdcFullScaleStack");
  else 
    theStackADC_=255;
  if ( conf_.exists("FirstStackLayer") ) theFirstStack_=conf_.getParameter<int>("FirstStackLayer");
  else
    theFirstStack_=5;
  
  // Get the constants for the miss-calibration studies
  doMissCalibrate=conf_.getUntrackedParameter<bool>("MissCalibrate",true); 
  doSplitClusters = conf.getParameter<bool>("SplitClusters");
  theBuffer.setSize( theNumOfRows, theNumOfCols );
}

PixelThresholdClusterizer::~PixelThresholdClusterizer ( )

Definition at line 71 of file PixelThresholdClusterizer.cc.

{}

Member Function Documentation

int PixelThresholdClusterizer::calibrate	(	int	adc,
		int	col,
		int	row
	)		`[private]`

Definition at line 209 of file PixelThresholdClusterizer.cc.

References HI_PhotonSkim_cff::electrons, and PXBDetId::layer().

{
  int electrons = 0;
  int layer= 0;
  if (DetId(detid_).subdetId()==1){ layer = PXBDetId(detid_).layer();}

  if ( doMissCalibrate ) 
    {
      // do not perform calibration if pixel is dead!
      
      if ( !theSiPixelGainCalibrationService_->isDead(detid_,col,row) && 
           !theSiPixelGainCalibrationService_->isNoisy(detid_,col,row) )
        {
          
          // Linear approximation of the TANH response
          // Pixel(0,0,0)
          //const float gain = 2.95; // 1 ADC = 2.95 VCALs (1/0.339)
          //const float pedestal = -83.; // -28/0.339
          // Roc-0 average
          //const float gain = 1./0.357; // 1 ADC = 2.80 VCALs 
          //const float pedestal = -28.2 * gain; // -79.
          
          float DBgain     = theSiPixelGainCalibrationService_->getGain(detid_, col, row);
          float DBpedestal = theSiPixelGainCalibrationService_->getPedestal(detid_, col, row) * DBgain;
          
          
          // Roc-6 average
          //const float gain = 1./0.313; // 1 ADC = 3.19 VCALs 
          //const float pedestal = -6.2 * gain; // -19.8
          // 
          float vcal = adc * DBgain - DBpedestal;
          
          // atanh calibration 
          // Roc-6 average
          //const float p0 = 0.00492;
          //const float p1 = 1.998;
          //const float p2 = 90.6;
          //const float p3 = 134.1; 
          // Roc-6 average
          //const float p0 = 0.00382;
          //const float p1 = 0.886;
          //const float p2 = 112.7;
          //const float p3 = 113.0; 
          //float vcal = ( atanh( (adc-p3)/p2) + p1)/p0;
          
          electrons = int( vcal * theConversionFactor + theOffset); 
        }
    }
  else 
    { // No misscalibration in the digitizer
      // Simple (default) linear gain 
      const float gain = 135.; // 1 ADC = 135 electrons
      const float pedestal = 0.; //
      electrons = int(adc * gain + pedestal);
      if (layer>=theFirstStack_) {
        if (theStackADC_==1&&adc==1)
          {
            electrons = int(255*135); // Arbitrarily use overflow value.
          }
        if (theStackADC_>1&&theStackADC_!=255&&adc>=1)
          {
            const float gain = 135.; // 1 ADC = 135 electrons
            electrons = int((adc-1) * gain * 255/float(theStackADC_-1));
          }
      }
    }
  
  return electrons;
}

void PixelThresholdClusterizer::clear_buffer	(	DigiIterator	begin,
		DigiIterator	end
	)		`[private]`

Clear the internal buffer array.

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

TO DO: ask Danek... wouldn't it be faster to simply memcopy() zeros into the whole buffer array?

Definition at line 177 of file PixelThresholdClusterizer.cc.

References end.

{
  for(DigiIterator di = begin; di != end; ++di ) 
    {
      theBuffer.set_adc( di->row(), di->column(), 0 );   // reset pixel adc to 0
    }
}

void PixelThresholdClusterizer::clusterizeDetUnit	(	const edm::DetSet< PixelDigi > &	input,
		const PixelGeomDetUnit *	pixDet,
		const std::vector< short > &	badChannels,
		edmNew::DetSetVector< SiPixelCluster >::FastFiller &	output
	)		`[virtual]`

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

Implements PixelClusterizerBase.

Definition at line 111 of file PixelThresholdClusterizer.cc.

References edm::DetSet< T >::begin(), begin, SiPixelCluster::charge(), edm::DetSet< T >::detId(), edm::DetSet< T >::end(), end, i, edmNew::DetSetVector< T >::FastFiller::push_back(), and HcalObjRepresent::setup().

                                                                                                         {
  
  DigiIterator begin = input.begin();
  DigiIterator end   = input.end();
  
  // Do not bother for empty detectors
  //if (begin == end) cout << " PixelThresholdClusterizer::clusterizeDetUnit - No digis to clusterize";
  
  //  Set up the clusterization on this DetId.
  if ( !setup(pixDet) ) 
    return;
  
  detid_ = input.detId();
  
  //  Copy PixelDigis to the buffer array; select the seed pixels
  //  on the way, and store them in theSeeds.
  copy_to_buffer(begin, end);
  
  //  At this point we know the number of seeds on this DetUnit, and thus
  //  also the maximal number of possible clusters, so resize theClusters
  //  in order to make vector<>::push_back() efficient.
  // output.reserve ( theSeeds.size() ); //GPetruc: It is better *not* to reserve, with the new DetSetVector!
  
  
  //  Loop over all seeds.  TO DO: wouldn't using iterators be faster?
  //  edm::LogError("PixelThresholdClusterizer") <<  "Starting clusterizing" << endl;
  for (unsigned int i = 0; i < theSeeds.size(); i++) 
    {
      
      // Gavril : The charge of seeds that were already inlcuded in clusters is set to 1 electron
      // so we don't want to call "make_cluster" for these cases 
      if ( theBuffer(theSeeds[i]) >= theSeedThreshold ) 
        {  // Is this seed still valid?
          //  Make a cluster around this seed
          SiPixelCluster cluster = make_cluster( theSeeds[i] , output);
          
          //  Check if the cluster is above threshold  
          // (TO DO: one is signed, other unsigned, gcc warns...)
          if ( cluster.charge() >= theClusterThreshold) 
            {
              //        cout << "putting in this cluster" << endl;
              output.push_back( cluster );
            }
        }
    }
  
  // Erase the seeds.
  theSeeds.clear();
  
  //  Need to clean unused pixels from the buffer array.
  clear_buffer(begin, end);
  
}

void PixelThresholdClusterizer::copy_to_buffer	(	DigiIterator	begin,
		DigiIterator	end
	)		`[private]`

Copy adc counts from PixelDigis into the buffer, identify seeds.

Definition at line 188 of file PixelThresholdClusterizer.cc.

References ecalMGPA::adc(), and end.

{
  for(DigiIterator di = begin; di != end; ++di) 
    {
      int row = di->row();
      int col = di->column();
      int adc = calibrate(di->adc(),col,row); // convert ADC -> electrons
      if ( adc >= thePixelThreshold) 
        {
          theBuffer.set_adc( row, col, adc);
          if ( adc >= theSeedThreshold) 
            { 
              theSeeds.push_back( SiPixelCluster::PixelPos(row,col) );
            }
        }
    }
}

SiPixelCluster PixelThresholdClusterizer::make_cluster	(	const SiPixelCluster::PixelPos &	pix,
		edmNew::DetSetVector< SiPixelCluster >::FastFiller &	output
	)		`[private]`

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

Definition at line 315 of file PixelThresholdClusterizer.cc.

References SiPixelCluster::add(), trackerHits::c, SiPixelCluster::charge(), SiPixelCluster::PixelPos::col(), i, SiPixelCluster::pixels(), edmNew::DetSetVector< T >::FastFiller::push_back(), alignCSCRings::r, and SiPixelCluster::PixelPos::row().

{
  
  //First we acquire the seeds for the clusters
  int seed_adc;
  stack<SiPixelCluster::PixelPos, vector<SiPixelCluster::PixelPos> > dead_pixel_stack;
  
  //The individual modules have been loaded into a buffer.
  //After each pixel has been considered by the clusterizer, we set the adc count to 1
  //to mark that we have already considered it.
  //The only difference between dead/noisy pixels and standard ones is that for dead/noisy pixels,
  //We consider the charge of the pixel to always be zero.

  if ( doMissCalibrate &&
       (theSiPixelGainCalibrationService_->isDead(detid_,pix.col(),pix.row()) || 
        theSiPixelGainCalibrationService_->isNoisy(detid_,pix.col(),pix.row())) )
    {
      seed_adc = 0;
      theBuffer.set_adc(pix, 1);
    }
  else
    {
      seed_adc = theBuffer(pix.row(), pix.col());
      theBuffer.set_adc( pix, 1);
    }
  
  AccretionCluster acluster;
  acluster.add(pix, seed_adc);
  
  //Here we search all pixels adjacent to all pixels in the cluster.
  bool dead_flag = false;
  while ( ! acluster.empty()) 
    {
      //This is the standard algorithm to find and add a pixel
      auto curInd = acluster.top(); acluster.pop();
      for ( auto r = acluster.x[curInd]-1; r <= acluster.x[curInd]+1; ++r) 
        {
          for ( auto c = acluster.y[curInd]-1; c <= acluster.y[curInd]+1; ++c) 
            {
              if ( theBuffer(r,c) >= thePixelThreshold) 
                {
                  
                  SiPixelCluster::PixelPos newpix(r,c);
                  if (!acluster.add( newpix, theBuffer(r,c))) goto endClus;
                  theBuffer.set_adc( newpix, 1);
                }
             

              /* //Commenting out the addition of dead pixels to the cluster until further testing -- dfehling 06/09
              //Check on the bounds of the module; this is to keep the isDead and isNoisy modules from returning errors 
              else if(r>= 0 && c >= 0 && (r <= (theNumOfRows-1.)) && (c <= (theNumOfCols-1.))){ 
              //Check for dead/noisy pixels check that the buffer is not -1 (already considered).  Check whether we want to split clusters separated by dead pixels or not.
              if((theSiPixelGainCalibrationService_->isDead(detid_,c,r) || theSiPixelGainCalibrationService_->isNoisy(detid_,c,r)) && theBuffer(r,c) != 1){
              
              //If a pixel is dead or noisy, check to see if we want to split the clusters or not.  
              //Push it into a dead pixel stack in case we want to split the clusters.  Otherwise add it to the cluster.
              //If we are splitting the clusters, we will iterate over the dead pixel stack later.
              
              SiPixelCluster::PixelPos newpix(r,c);
              if(!doSplitClusters){
              
              cluster.add(newpix, theBuffer(r,c));}
              else if(doSplitClusters){
              dead_pixel_stack.push(newpix);
              dead_flag = true;}
              
              theBuffer.set_adc(newpix, 1);
              } 
              
              }
              */
            


            }
        }
      
    }  // while accretion
 endClus:
  SiPixelCluster cluster(acluster.isize,acluster.adc, acluster.x,acluster.y, acluster.xmin,acluster.ymin);
  //Here we split the cluster, if the flag to do so is set and we have found a dead or noisy pixel.
  
  if (dead_flag && doSplitClusters) 
    {
      //Set the first cluster equal to the existing cluster.
      SiPixelCluster first_cluster = cluster;
      bool have_second_cluster = false;
      while ( !dead_pixel_stack.empty() )
        {
          //consider each found dead pixel
          SiPixelCluster::PixelPos deadpix = dead_pixel_stack.top(); dead_pixel_stack.pop();
          theBuffer.set_adc(deadpix, 1);
         
          //Clusterize the split cluster using the dead pixel as a seed
          SiPixelCluster second_cluster = make_cluster(deadpix, output);
          
          //If both clusters would normally have been found by the clusterizer, put them into output
          if ( second_cluster.charge() >= theClusterThreshold && 
               first_cluster.charge() >= theClusterThreshold )
            {
              output.push_back( second_cluster );
              have_second_cluster = true;       
            }
          
          //We also want to keep the merged cluster in data and let the RecHit algorithm decide which set to keep
          //This loop adds the second cluster to the first.
          const std::vector<SiPixelCluster::Pixel>& branch_pixels = second_cluster.pixels();
          for ( unsigned int i = 0; i<branch_pixels.size(); i++)
            {
              int temp_x = branch_pixels[i].x;
              int temp_y = branch_pixels[i].y;
              int temp_adc = branch_pixels[i].adc;
              SiPixelCluster::PixelPos newpix(temp_x, temp_y);
              cluster.add(newpix, temp_adc);}
        }
      
      //Remember to also add the first cluster if we added the second one.
      if ( first_cluster.charge() >= theClusterThreshold && have_second_cluster) 
        {
          output.push_back( first_cluster );
        }
    }
  
  return cluster;
}

bool PixelThresholdClusterizer::setup ( const PixelGeomDetUnit * pixDet ) [private]

Private helper methods:

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 77 of file PixelThresholdClusterizer.cc.

References PixelTopology::ncolumns(), PixelTopology::nrows(), and PixelGeomDetUnit::specificTopology().

{
  // Cache the topology.
  const PixelTopology & topol = pixDet->specificTopology();
  
  // Get the new sizes.
  int nrows = topol.nrows();      // rows in x
  int ncols = topol.ncolumns();   // cols in y
  
  theNumOfRows = nrows;  // Set new sizes
  theNumOfCols = ncols;
  
  if ( nrows > theBuffer.rows() || 
       ncols > theBuffer.columns() ) 
    { // change only when a larger is needed
      //if( nrows != theNumOfRows || ncols != theNumOfCols ) {
      //cout << " PixelThresholdClusterizer: pixel buffer redefined to " 
      // << nrows << " * " << ncols << endl;      
      //theNumOfRows = nrows;  // Set new sizes
      //theNumOfCols = ncols;
      // Resize the buffer
      theBuffer.setSize(nrows,ncols);  // Modify
      bufferAlreadySet = true;
    }
  
  return true;   
}