Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-05-31 22:26:01

0001 //----------------------------------------------------------------------------
0002 //! \class PixelThresholdClusterizer
0003 //! \brief A specific threshold-based pixel clustering algorithm
0004 //!
0005 //! General logic of PixelThresholdClusterizer:
0006 //!
0007 //! The clusterization is performed on a matrix with size
0008 //! equal to the size of the pixel detector, each cell containing
0009 //! the ADC count of the corresponding pixel.
0010 //! The matrix is reset after each clusterization.
0011 //!
0012 //! The search starts from seed pixels, i.e. pixels with sufficiently
0013 //! large amplitudes, found at the time of filling of the matrix
0014 //! and stored in a SiPixelArrayBuffer.
0015 //!
0016 //! Translate the pixel charge to electrons, we are suppose to
0017 //! do the calibrations ADC->electrons here.
0018 //! Modify the thresholds to be in electrons, convert adc to electrons. d.k. 20/3/06
0019 //! Get rid of the noiseVector. d.k. 28/3/06
0020 //----------------------------------------------------------------------------
0021 
0022 // Our own includes
0023 #include "PixelThresholdClusterizer.h"
0024 #include "SiPixelArrayBuffer.h"
0025 #include "CondFormats/SiPixelObjects/interface/SiPixelGainCalibrationOffline.h"
0026 // Geometry
0027 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0028 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0029 //#include "Geometry/CommonTopologies/RectangularPixelTopology.h"
0030 
0031 // STL
0032 #include <stack>
0033 #include <vector>
0034 #include <iostream>
0035 #include <atomic>
0036 #include <algorithm>
0037 #include <limits>
0038 
0039 //----------------------------------------------------------------------------
0040 //! Constructor:
0041 //!  Initilize the buffer to hold pixels from a detector module.
0042 //!  This is a vector of 44k ints, stays valid all the time.
0043 //----------------------------------------------------------------------------
0044 PixelThresholdClusterizer::PixelThresholdClusterizer(edm::ParameterSet const& conf)
0045     :  // Get thresholds in electrons
0046       thePixelThreshold(conf.getParameter<int>("ChannelThreshold")),
0047       theSeedThreshold(conf.getParameter<int>("SeedThreshold")),
0048       theClusterThreshold(conf.getParameter<int>("ClusterThreshold")),
0049       theClusterThreshold_L1(conf.getParameter<int>("ClusterThreshold_L1")),
0050       theConversionFactor(conf.getParameter<int>("VCaltoElectronGain")),
0051       theConversionFactor_L1(conf.getParameter<int>("VCaltoElectronGain_L1")),
0052       theOffset(conf.getParameter<int>("VCaltoElectronOffset")),
0053       theOffset_L1(conf.getParameter<int>("VCaltoElectronOffset_L1")),
0054       theElectronPerADCGain(conf.getParameter<double>("ElectronPerADCGain")),
0055       doPhase2Calibration(conf.getParameter<bool>("Phase2Calibration")),
0056       thePhase2ReadoutMode(conf.getParameter<int>("Phase2ReadoutMode")),
0057       thePhase2DigiBaseline(conf.getParameter<double>("Phase2DigiBaseline")),
0058       thePhase2KinkADC(conf.getParameter<int>("Phase2KinkADC")),
0059       theNumOfRows(0),
0060       theNumOfCols(0),
0061       theDetid(0),
0062       // Get the constants for the miss-calibration studies
0063       doMissCalibrate(conf.getParameter<bool>("MissCalibrate")),
0064       doSplitClusters(conf.getParameter<bool>("SplitClusters")) {
0065   theBuffer.setSize(theNumOfRows, theNumOfCols);
0066   theFakePixels.clear();
0067 }
0068 /////////////////////////////////////////////////////////////////////////////
0069 PixelThresholdClusterizer::~PixelThresholdClusterizer() {}
0070 
0071 // Configuration descriptions
0072 void PixelThresholdClusterizer::fillPSetDescription(edm::ParameterSetDescription& desc) {
0073   desc.add<int>("ChannelThreshold", 1000);
0074   desc.add<bool>("MissCalibrate", true);
0075   desc.add<bool>("SplitClusters", false);
0076   desc.add<int>("VCaltoElectronGain", 65);
0077   desc.add<int>("VCaltoElectronGain_L1", 65);
0078   desc.add<int>("VCaltoElectronOffset", -414);
0079   desc.add<int>("VCaltoElectronOffset_L1", -414);
0080   desc.add<int>("SeedThreshold", 1000);
0081   desc.add<int>("ClusterThreshold_L1", 4000);
0082   desc.add<int>("ClusterThreshold", 4000);
0083   desc.add<double>("ElectronPerADCGain", 135.);
0084   desc.add<bool>("Phase2Calibration", false);
0085   desc.add<int>("Phase2ReadoutMode", -1);
0086   desc.add<double>("Phase2DigiBaseline", 1200.);
0087   desc.add<int>("Phase2KinkADC", 8);
0088 }
0089 
0090 //----------------------------------------------------------------------------
0091 //!  Prepare the Clusterizer to work on a particular DetUnit.  Re-init the
0092 //!  size of the panel/plaquette (so update nrows and ncols),
0093 //----------------------------------------------------------------------------
0094 bool PixelThresholdClusterizer::setup(const PixelGeomDetUnit* pixDet) {
0095   // Cache the topology.
0096   const PixelTopology& topol = pixDet->specificTopology();
0097 
0098   // Get the new sizes.
0099   int nrows = topol.nrows();     // rows in x
0100   int ncols = topol.ncolumns();  // cols in y
0101 
0102   theNumOfRows = nrows;  // Set new sizes
0103   theNumOfCols = ncols;
0104 
0105   if (nrows > theBuffer.rows() || ncols > theBuffer.columns()) {  // change only when a larger is needed
0106     if (nrows != theNumOfRows || ncols != theNumOfCols)
0107       edm::LogWarning("setup()") << "pixel buffer redefined to" << nrows << " * " << ncols;
0108     //theNumOfRows = nrows;  // Set new sizes
0109     //theNumOfCols = ncols;
0110     // Resize the buffer
0111     theBuffer.setSize(nrows, ncols);  // Modify
0112   }
0113 
0114   theFakePixels.resize(nrows * ncols, false);
0115 
0116   return true;
0117 }
0118 //----------------------------------------------------------------------------
0119 //!  \brief Cluster pixels.
0120 //!  This method operates on a matrix of pixels
0121 //!  and finds the largest contiguous cluster around
0122 //!  each seed pixel.
0123 //!  Input and output data stored in DetSet
0124 //----------------------------------------------------------------------------
0125 template <typename T>
0126 void PixelThresholdClusterizer::clusterizeDetUnitT(const T& input,
0127                                                    const PixelGeomDetUnit* pixDet,
0128                                                    const TrackerTopology* tTopo,
0129                                                    const std::vector<short>& badChannels,
0130                                                    edmNew::DetSetVector<SiPixelCluster>::FastFiller& output) {
0131   typename T::const_iterator begin = input.begin();
0132   typename T::const_iterator end = input.end();
0133 
0134   // this should never happen and the raw2digi does not create empty detsets
0135   if (begin == end) {
0136     edm::LogError("PixelThresholdClusterizer") << "@SUB=PixelThresholdClusterizer::clusterizeDetUnitT()"
0137                                                << " No digis to clusterize";
0138   }
0139 
0140   //  Set up the clusterization on this DetId.
0141   if (!setup(pixDet))
0142     return;
0143 
0144   theDetid = input.detId();
0145 
0146   // Set separate cluster threshold for L1 (needed for phase1)
0147   auto clusterThreshold = theClusterThreshold;
0148   theLayer = (DetId(theDetid).subdetId() == 1) ? tTopo->pxbLayer(theDetid) : 0;
0149   if (theLayer == 1)
0150     clusterThreshold = theClusterThreshold_L1;
0151 
0152   //  Copy PixelDigis to the buffer array; select the seed pixels
0153   //  on the way, and store them in theSeeds.
0154   if (end > begin)
0155     copy_to_buffer(begin, end);
0156 
0157   assert(output.empty());
0158   //  Loop over all seeds.  TO DO: wouldn't using iterators be faster?
0159   for (unsigned int i = 0; i < theSeeds.size(); i++) {
0160     // Gavril : The charge of seeds that were already inlcuded in clusters is set to 1 electron
0161     // so we don't want to call "make_cluster" for these cases
0162     if (theBuffer(theSeeds[i]) >= theSeedThreshold) {  // Is this seed still valid?
0163       //  Make a cluster around this seed
0164       SiPixelCluster&& cluster = make_cluster(theSeeds[i], output);
0165 
0166       //  Check if the cluster is above threshold
0167       // (TO DO: one is signed, other unsigned, gcc warns...)
0168       if (cluster.charge() >= clusterThreshold) {
0169         // sort by row (x)
0170         output.push_back(std::move(cluster));
0171         std::push_heap(output.begin(), output.end(), [](SiPixelCluster const& cl1, SiPixelCluster const& cl2) {
0172           return cl1.minPixelRow() < cl2.minPixelRow();
0173         });
0174       }
0175     }
0176   }
0177   // sort by row (x)   maybe sorting the seed would suffice....
0178   std::sort_heap(output.begin(), output.end(), [](SiPixelCluster const& cl1, SiPixelCluster const& cl2) {
0179     return cl1.minPixelRow() < cl2.minPixelRow();
0180   });
0181 
0182   // Erase the seeds.
0183   theSeeds.clear();
0184 
0185   //  Need to clean unused pixels from the buffer array.
0186   clear_buffer(begin, end);
0187 
0188   theFakePixels.clear();
0189 }
0190 
0191 //----------------------------------------------------------------------------
0192 //!  \brief Clear the internal buffer array.
0193 //!
0194 //!  Pixels which are not part of recognized clusters are NOT ERASED
0195 //!  during the cluster finding.  Erase them now.
0196 //!
0197 //!  TO DO: ask Danek... wouldn't it be faster to simply memcopy() zeros into
0198 //!  the whole buffer array?
0199 //----------------------------------------------------------------------------
0200 void PixelThresholdClusterizer::clear_buffer(DigiIterator begin, DigiIterator end) {
0201   for (DigiIterator di = begin; di != end; ++di) {
0202     theBuffer.set_adc(di->row(), di->column(), 0);  // reset pixel adc to 0
0203   }
0204 }
0205 
0206 void PixelThresholdClusterizer::clear_buffer(ClusterIterator begin, ClusterIterator end) {
0207   for (ClusterIterator ci = begin; ci != end; ++ci) {
0208     for (int i = 0; i < ci->size(); ++i) {
0209       const SiPixelCluster::Pixel pixel = ci->pixel(i);
0210 
0211       theBuffer.set_adc(pixel.x, pixel.y, 0);  // reset pixel adc to 0
0212     }
0213   }
0214 }
0215 
0216 //----------------------------------------------------------------------------
0217 //! \brief Copy adc counts from PixelDigis into the buffer, identify seeds.
0218 //----------------------------------------------------------------------------
0219 void PixelThresholdClusterizer::copy_to_buffer(DigiIterator begin, DigiIterator end) {
0220 #ifdef PIXELREGRESSION
0221   static std::atomic<int> s_ic = 0;
0222   in ic = ++s_ic;
0223   if (ic == 1) {
0224     // std::cout << (doMissCalibrate ? "VI from db" : "VI linear") << std::endl;
0225   }
0226 #endif
0227 
0228   //If called with empty/invalid DetSet, warn the user
0229   if (end <= begin) {
0230     edm::LogWarning("PixelThresholdClusterizer") << " copy_to_buffer called with empty or invalid range" << std::endl;
0231     return;
0232   }
0233 
0234   int electron[end - begin];  // pixel charge in electrons
0235   memset(electron, 0, (end - begin) * sizeof(int));
0236 
0237   if (doPhase2Calibration) {
0238     int i = 0;
0239     for (DigiIterator di = begin; di != end; ++di) {
0240       electron[i] = calibrate(di->adc(), di->column(), di->row());
0241       i++;
0242     }
0243     assert(i == (end - begin));
0244   }
0245 
0246   else {
0247     if (doMissCalibrate) {
0248       if (theLayer == 1) {
0249         (*theSiPixelGainCalibrationService_)
0250             .calibrate(theDetid, begin, end, theConversionFactor_L1, theOffset_L1, electron);
0251       } else {
0252         (*theSiPixelGainCalibrationService_).calibrate(theDetid, begin, end, theConversionFactor, theOffset, electron);
0253       }
0254     } else {
0255       int i = 0;
0256       const float gain = theElectronPerADCGain;  // default: 1 ADC = 135 electrons
0257       for (DigiIterator di = begin; di != end; ++di) {
0258         auto adc = di->adc();
0259         const float pedestal = 0.;  //
0260         electron[i] = int(adc * gain + pedestal);
0261         ++i;
0262       }
0263       assert(i == (end - begin));
0264     }
0265   }
0266 
0267   int i = 0;
0268 #ifdef PIXELREGRESSION
0269   static std::atomic<int> eqD = 0;
0270 #endif
0271   for (DigiIterator di = begin; di != end; ++di) {
0272     int row = di->row();
0273     int col = di->column();
0274     // VV: do not calibrate a fake pixel, it already has a unit of 10e-:
0275     int adc = (di->flag() != 0) ? di->adc() * 10 : electron[i];  // this is in electrons
0276     i++;
0277 
0278 #ifdef PIXELREGRESSION
0279     int adcOld = calibrate(di->adc(), col, row);
0280     //assert(adc==adcOld);
0281     if (adc != adcOld)
0282       std::cout << "VI " << eqD << ' ' << ic << ' ' << end - begin << ' ' << i << ' ' << di->adc() << ' ' << adc << ' '
0283                 << adcOld << std::endl;
0284     else
0285       ++eqD;
0286 #endif
0287 
0288     if (adc < 100)
0289       adc = 100;  // put all negative pixel charges into the 100 elec bin
0290     /* This is semi-random good number. The exact number (in place of 100) is irrelevant from the point 
0291        of view of the final cluster charge since these are typically >= 20000.
0292     */
0293 
0294     if (adc >= thePixelThreshold) {
0295       theBuffer.set_adc(row, col, adc);
0296       // VV: add pixel to the fake list. Only when running on digi collection
0297       if (di->flag() != 0)
0298         theFakePixels[row * theNumOfCols + col] = true;
0299       if (adc >= theSeedThreshold)
0300         theSeeds.push_back(SiPixelCluster::PixelPos(row, col));
0301     }
0302   }
0303   assert(i == (end - begin));
0304 }
0305 
0306 void PixelThresholdClusterizer::copy_to_buffer(ClusterIterator begin, ClusterIterator end) {
0307   // loop over clusters
0308   for (ClusterIterator ci = begin; ci != end; ++ci) {
0309     // loop over pixels
0310     for (int i = 0; i < ci->size(); ++i) {
0311       const SiPixelCluster::Pixel pixel = ci->pixel(i);
0312 
0313       int row = pixel.x;
0314       int col = pixel.y;
0315       int adc = pixel.adc;
0316       if (adc >= thePixelThreshold) {
0317         theBuffer.add_adc(row, col, adc);
0318         if (adc >= theSeedThreshold)
0319           theSeeds.push_back(SiPixelCluster::PixelPos(row, col));
0320       }
0321     }
0322   }
0323 }
0324 
0325 //----------------------------------------------------------------------------
0326 // Calibrate adc counts to electrons
0327 //-----------------------------------------------------------------
0328 int PixelThresholdClusterizer::calibrate(int adc, int col, int row) {
0329   int electrons = 0;
0330 
0331   if (doPhase2Calibration) {
0332     const float gain = theElectronPerADCGain;
0333     int p2rm = (thePhase2ReadoutMode < -1 ? -1 : thePhase2ReadoutMode);
0334 
0335     if (p2rm == -1) {
0336       electrons = int(adc * gain);
0337     } else {
0338       if (adc < thePhase2KinkADC) {
0339         electrons = int((adc + 0.5) * gain);
0340       } else {
0341         const int dualslopeparam = (thePhase2ReadoutMode < 10 ? thePhase2ReadoutMode : 10);
0342         const int dualslope = int(dualslopeparam <= 1 ? 1. : pow(2, dualslopeparam - 1));
0343         adc -= thePhase2KinkADC;
0344         adc *= dualslope;
0345         adc += thePhase2KinkADC;
0346         electrons = int((adc + 0.5 * dualslope) * gain);
0347       }
0348       electrons += int(thePhase2DigiBaseline);
0349     }
0350 
0351     return electrons;
0352   }
0353 
0354   if (doMissCalibrate) {
0355     // do not perform calibration if pixel is dead!
0356 
0357     if (!theSiPixelGainCalibrationService_->isDead(theDetid, col, row) &&
0358         !theSiPixelGainCalibrationService_->isNoisy(theDetid, col, row)) {
0359       // Linear approximation of the TANH response
0360       // Pixel(0,0,0)
0361       //const float gain = 2.95; // 1 ADC = 2.95 VCALs (1/0.339)
0362       //const float pedestal = -83.; // -28/0.339
0363       // Roc-0 average
0364       //const float gain = 1./0.357; // 1 ADC = 2.80 VCALs
0365       //const float pedestal = -28.2 * gain; // -79.
0366 
0367       float DBgain = theSiPixelGainCalibrationService_->getGain(theDetid, col, row);
0368       float pedestal = theSiPixelGainCalibrationService_->getPedestal(theDetid, col, row);
0369       float DBpedestal = pedestal * DBgain;
0370 
0371       // Roc-6 average
0372       //const float gain = 1./0.313; // 1 ADC = 3.19 VCALs
0373       //const float pedestal = -6.2 * gain; // -19.8
0374       //
0375       float vcal = adc * DBgain - DBpedestal;
0376 
0377       // atanh calibration
0378       // Roc-6 average
0379       //const float p0 = 0.00492;
0380       //const float p1 = 1.998;
0381       //const float p2 = 90.6;
0382       //const float p3 = 134.1;
0383       // Roc-6 average
0384       //const float p0 = 0.00382;
0385       //const float p1 = 0.886;
0386       //const float p2 = 112.7;
0387       //const float p3 = 113.0;
0388       //float vcal = ( atanh( (adc-p3)/p2) + p1)/p0;
0389 
0390       if (theLayer == 1) {
0391         electrons = int(vcal * theConversionFactor_L1 + theOffset_L1);
0392       } else {
0393         electrons = int(vcal * theConversionFactor + theOffset);
0394       }
0395     }
0396   } else {  // No misscalibration in the digitizer
0397     // Simple (default) linear gain
0398     const float gain = theElectronPerADCGain;  // default: 1 ADC = 135 electrons
0399     const float pedestal = 0.;                 //
0400     electrons = int(adc * gain + pedestal);
0401   }
0402 
0403   return electrons;
0404 }
0405 
0406 //----------------------------------------------------------------------------
0407 //!  \brief The actual clustering algorithm: group the neighboring pixels around the seed.
0408 //----------------------------------------------------------------------------
0409 SiPixelCluster PixelThresholdClusterizer::make_cluster(const SiPixelCluster::PixelPos& pix,
0410                                                        edmNew::DetSetVector<SiPixelCluster>::FastFiller& output) {
0411   //First we acquire the seeds for the clusters
0412   uint16_t seed_adc;
0413   std::stack<SiPixelCluster::PixelPos, std::vector<SiPixelCluster::PixelPos> > dead_pixel_stack;
0414 
0415   //The individual modules have been loaded into a buffer.
0416   //After each pixel has been considered by the clusterizer, we set the adc count to 1
0417   //to mark that we have already considered it.
0418   //The only difference between dead/noisy pixels and standard ones is that for dead/noisy pixels,
0419   //We consider the charge of the pixel to always be zero.
0420 
0421   /*  this is not possible as dead and noisy pixel cannot make it into a seed...
0422   if ( doMissCalibrate &&
0423        (theSiPixelGainCalibrationService_->isDead(theDetid,pix.col(),pix.row()) || 
0424     theSiPixelGainCalibrationService_->isNoisy(theDetid,pix.col(),pix.row())) )
0425     {
0426       std::cout << "IMPOSSIBLE" << std::endl;
0427       seed_adc = 0;
0428       theBuffer.set_adc(pix, 1);
0429     }
0430     else {
0431   */
0432   // Note: each ADC value is limited here to 65535 (std::numeric_limits<uint16_t>::max),
0433   //       as it is later stored as uint16_t in SiPixelCluster and PixelClusterizerBase/AccretionCluster
0434   //       (reminder: ADC values here may be expressed in number of electrons)
0435   seed_adc = std::min(theBuffer(pix.row(), pix.col()), int(std::numeric_limits<uint16_t>::max()));
0436   theBuffer.set_adc(pix, 1);
0437   //  }
0438 
0439   AccretionCluster acluster, cldata;
0440   acluster.add(pix, seed_adc);
0441   cldata.add(pix, seed_adc);
0442 
0443   //Here we search all pixels adjacent to all pixels in the cluster.
0444   bool dead_flag = false;
0445   while (!acluster.empty()) {
0446     //This is the standard algorithm to find and add a pixel
0447     auto curInd = acluster.top();
0448     acluster.pop();
0449     for (auto c = std::max(0, int(acluster.y[curInd]) - 1);
0450          c < std::min(int(acluster.y[curInd]) + 2, theBuffer.columns());
0451          ++c) {
0452       for (auto r = std::max(0, int(acluster.x[curInd]) - 1);
0453            r < std::min(int(acluster.x[curInd]) + 2, theBuffer.rows());
0454            ++r) {
0455         if (theBuffer(r, c) >= thePixelThreshold) {
0456           SiPixelCluster::PixelPos newpix(r, c);
0457           auto const newpix_adc = std::min(theBuffer(r, c), int(std::numeric_limits<uint16_t>::max()));
0458           if (!acluster.add(newpix, newpix_adc))
0459             goto endClus;
0460           // VV: no fake pixels in cluster, leads to non-contiguous clusters
0461           if (!theFakePixels[r * theNumOfCols + c]) {
0462             cldata.add(newpix, newpix_adc);
0463           }
0464           theBuffer.set_adc(newpix, 1);
0465         }
0466 
0467         /* //Commenting out the addition of dead pixels to the cluster until further testing -- dfehling 06/09
0468           //Check on the bounds of the module; this is to keep the isDead and isNoisy modules from returning errors 
0469           else if(r>= 0 && c >= 0 && (r <= (theNumOfRows-1.)) && (c <= (theNumOfCols-1.))){ 
0470           //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.
0471           if((theSiPixelGainCalibrationService_->isDead(theDetid,c,r) || theSiPixelGainCalibrationService_->isNoisy(theDetid,c,r)) && theBuffer(r,c) != 1){
0472           
0473           //If a pixel is dead or noisy, check to see if we want to split the clusters or not.  
0474           //Push it into a dead pixel stack in case we want to split the clusters.  Otherwise add it to the cluster.
0475           //If we are splitting the clusters, we will iterate over the dead pixel stack later.
0476           
0477           SiPixelCluster::PixelPos newpix(r,c);
0478           if(!doSplitClusters){
0479 
0480           cluster.add(newpix, std::min(theBuffer(r, c), int(std::numeric_limits<uint16_t>::max())));}
0481           else if(doSplitClusters){
0482           dead_pixel_stack.push(newpix);
0483           dead_flag = true;}
0484           
0485           theBuffer.set_adc(newpix, 1);
0486           } 
0487           
0488           }
0489           */
0490       }
0491     }
0492 
0493   }  // while accretion
0494 endClus:
0495   SiPixelCluster cluster(cldata.isize, cldata.adc, cldata.x, cldata.y, cldata.xmin, cldata.ymin);
0496   //Here we split the cluster, if the flag to do so is set and we have found a dead or noisy pixel.
0497 
0498   if (dead_flag && doSplitClusters) {
0499     // Set separate cluster threshold for L1 (needed for phase1)
0500     auto clusterThreshold = theClusterThreshold;
0501     if (theLayer == 1)
0502       clusterThreshold = theClusterThreshold_L1;
0503 
0504     //Set the first cluster equal to the existing cluster.
0505     SiPixelCluster first_cluster = cluster;
0506     bool have_second_cluster = false;
0507     while (!dead_pixel_stack.empty()) {
0508       //consider each found dead pixel
0509       SiPixelCluster::PixelPos deadpix = dead_pixel_stack.top();
0510       dead_pixel_stack.pop();
0511       theBuffer.set_adc(deadpix, 1);
0512 
0513       //Clusterize the split cluster using the dead pixel as a seed
0514       SiPixelCluster second_cluster = make_cluster(deadpix, output);
0515 
0516       //If both clusters would normally have been found by the clusterizer, put them into output
0517       if (second_cluster.charge() >= clusterThreshold && first_cluster.charge() >= clusterThreshold) {
0518         output.push_back(second_cluster);
0519         have_second_cluster = true;
0520       }
0521 
0522       //We also want to keep the merged cluster in data and let the RecHit algorithm decide which set to keep
0523       //This loop adds the second cluster to the first.
0524       const std::vector<SiPixelCluster::Pixel>& branch_pixels = second_cluster.pixels();
0525       for (unsigned int i = 0; i < branch_pixels.size(); i++) {
0526         auto const temp_x = branch_pixels[i].x;
0527         auto const temp_y = branch_pixels[i].y;
0528         auto const temp_adc = branch_pixels[i].adc;
0529         SiPixelCluster::PixelPos newpix(temp_x, temp_y);
0530         cluster.add(newpix, temp_adc);
0531       }
0532     }
0533 
0534     //Remember to also add the first cluster if we added the second one.
0535     if (first_cluster.charge() >= clusterThreshold && have_second_cluster) {
0536       output.push_back(first_cluster);
0537       std::push_heap(output.begin(), output.end(), [](SiPixelCluster const& cl1, SiPixelCluster const& cl2) {
0538         return cl1.minPixelRow() < cl2.minPixelRow();
0539       });
0540     }
0541   }
0542 
0543   return cluster;
0544 }