Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-03-17 11:19:34

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