File indexing completed on 2024-12-01 23:40:36
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023 #include "PixelThresholdClusterizer.h"
0024 #include "SiPixelArrayBuffer.h"
0025 #include "CondFormats/SiPixelObjects/interface/SiPixelGainCalibrationOffline.h"
0026
0027 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0028 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0029
0030
0031
0032 #include <stack>
0033 #include <vector>
0034 #include <iostream>
0035 #include <atomic>
0036 #include <algorithm>
0037 #include <limits>
0038
0039
0040
0041
0042
0043
0044 PixelThresholdClusterizer::PixelThresholdClusterizer(edm::ParameterSet const& conf)
0045 :
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
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
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
0095
0096
0097 bool PixelThresholdClusterizer::setup(const PixelGeomDetUnit* pixDet) {
0098
0099 const PixelTopology& topol = pixDet->specificTopology();
0100
0101
0102 int nrows = topol.nrows();
0103 int ncols = topol.ncolumns();
0104
0105 theNumOfRows = nrows;
0106 theNumOfCols = ncols;
0107
0108 if (nrows > theBuffer.rows() || ncols > theBuffer.columns()) {
0109 if (nrows != theNumOfRows || ncols != theNumOfCols)
0110 edm::LogWarning("setup()") << "pixel buffer redefined to" << nrows << " * " << ncols;
0111
0112
0113
0114 theBuffer.setSize(nrows, ncols);
0115 }
0116
0117 theFakePixels.resize(nrows * ncols, false);
0118
0119 thePixelOccurrence.resize(nrows * ncols, 0);
0120
0121 return true;
0122 }
0123
0124 #include "PixelThresholdClusterizer.icc"
0125
0126
0127
0128
0129
0130
0131
0132
0133
0134
0135 void PixelThresholdClusterizer::clear_buffer(DigiIterator begin, DigiIterator end) {
0136 for (DigiIterator di = begin; di != end; ++di) {
0137 theBuffer.set_adc(di->row(), di->column(), 0);
0138 }
0139 }
0140
0141 void PixelThresholdClusterizer::clear_buffer(ClusterIterator begin, ClusterIterator end) {
0142 for (ClusterIterator ci = begin; ci != end; ++ci) {
0143 for (int i = 0; i < ci->size(); ++i) {
0144 const SiPixelCluster::Pixel pixel = ci->pixel(i);
0145
0146 theBuffer.set_adc(pixel.x, pixel.y, 0);
0147 }
0148 }
0149 }
0150
0151
0152
0153
0154 void PixelThresholdClusterizer::copy_to_buffer(DigiIterator begin, DigiIterator end) {
0155 #ifdef PIXELREGRESSION
0156 static std::atomic<int> s_ic = 0;
0157 in ic = ++s_ic;
0158 if (ic == 1) {
0159
0160 }
0161 #endif
0162
0163
0164 if (end <= begin) {
0165 edm::LogWarning("PixelThresholdClusterizer") << " copy_to_buffer called with empty or invalid range" << std::endl;
0166 return;
0167 }
0168
0169 int electron[end - begin];
0170 memset(electron, 0, (end - begin) * sizeof(int));
0171
0172 if (doPhase2Calibration) {
0173 int i = 0;
0174 for (DigiIterator di = begin; di != end; ++di) {
0175 electron[i] = calibrate(di->adc(), di->column(), di->row());
0176 i++;
0177 }
0178 assert(i == (end - begin));
0179 }
0180
0181 else {
0182 if (doMissCalibrate) {
0183 if (theLayer == 1) {
0184 (*theSiPixelGainCalibrationService_)
0185 .calibrate(theDetid, begin, end, theConversionFactor_L1, theOffset_L1, electron);
0186 } else {
0187 (*theSiPixelGainCalibrationService_).calibrate(theDetid, begin, end, theConversionFactor, theOffset, electron);
0188 }
0189 } else {
0190 int i = 0;
0191 const float gain = theElectronPerADCGain;
0192 for (DigiIterator di = begin; di != end; ++di) {
0193 auto adc = di->adc();
0194 const float pedestal = 0.;
0195 electron[i] = int(adc * gain + pedestal);
0196 ++i;
0197 }
0198 assert(i == (end - begin));
0199 }
0200 }
0201
0202 int i = 0;
0203 #ifdef PIXELREGRESSION
0204 static std::atomic<int> eqD = 0;
0205 #endif
0206 for (DigiIterator di = begin; di != end; ++di) {
0207 int row = di->row();
0208 int col = di->column();
0209
0210 int adc = (di->flag() != 0) ? di->adc() * 10 : electron[i];
0211 i++;
0212
0213 #ifdef PIXELREGRESSION
0214 int adcOld = calibrate(di->adc(), col, row);
0215
0216 if (adc != adcOld)
0217 std::cout << "VI " << eqD << ' ' << ic << ' ' << end - begin << ' ' << i << ' ' << di->adc() << ' ' << adc << ' '
0218 << adcOld << std::endl;
0219 else
0220 ++eqD;
0221 #endif
0222
0223 if (adc < 100)
0224 adc = 100;
0225
0226
0227
0228
0229 thePixelOccurrence[theBuffer.index(row, col)]++;
0230 uint8_t occurrence =
0231 (!dropDuplicates) ? 1 : thePixelOccurrence[theBuffer.index(row, col)];
0232
0233 switch (occurrence) {
0234
0235 case 1:
0236 if (adc >= thePixelThreshold) {
0237 theBuffer.set_adc(row, col, adc);
0238
0239 if (di->flag() != 0)
0240 theFakePixels[row * theNumOfCols + col] = true;
0241 if (adc >= theSeedThreshold)
0242 theSeeds.push_back(SiPixelCluster::PixelPos(row, col));
0243 }
0244 break;
0245
0246
0247 case 2:
0248 theBuffer.set_adc(row, col, 0);
0249 std::remove(theSeeds.begin(), theSeeds.end(), SiPixelCluster::PixelPos(row, col));
0250 break;
0251
0252
0253 }
0254 }
0255 assert(i == (end - begin));
0256 }
0257
0258 void PixelThresholdClusterizer::copy_to_buffer(ClusterIterator begin, ClusterIterator end) {
0259
0260 for (ClusterIterator ci = begin; ci != end; ++ci) {
0261
0262 for (int i = 0; i < ci->size(); ++i) {
0263 const SiPixelCluster::Pixel pixel = ci->pixel(i);
0264
0265 int row = pixel.x;
0266 int col = pixel.y;
0267 int adc = pixel.adc;
0268 if (adc >= thePixelThreshold) {
0269 theBuffer.add_adc(row, col, adc);
0270 if (adc >= theSeedThreshold)
0271 theSeeds.push_back(SiPixelCluster::PixelPos(row, col));
0272 }
0273 }
0274 }
0275 }
0276
0277
0278
0279
0280 int PixelThresholdClusterizer::calibrate(int adc, int col, int row) {
0281 int electrons = 0;
0282
0283 if (doPhase2Calibration) {
0284 const float gain = theElectronPerADCGain;
0285 int p2rm = (thePhase2ReadoutMode < -1 ? -1 : thePhase2ReadoutMode);
0286
0287 if (p2rm == -1) {
0288 electrons = int(adc * gain);
0289 } else {
0290 if (adc < thePhase2KinkADC) {
0291 electrons = int((adc + 0.5) * gain);
0292 } else {
0293 const int dualslopeparam = (thePhase2ReadoutMode < 10 ? thePhase2ReadoutMode : 10);
0294 const int dualslope = int(dualslopeparam <= 1 ? 1. : pow(2, dualslopeparam - 1));
0295 adc -= thePhase2KinkADC;
0296 adc *= dualslope;
0297 adc += thePhase2KinkADC;
0298 electrons = int((adc + 0.5 * dualslope) * gain);
0299 }
0300 electrons += int(thePhase2DigiBaseline);
0301 }
0302
0303 return electrons;
0304 }
0305
0306 if (doMissCalibrate) {
0307
0308
0309 if (!theSiPixelGainCalibrationService_->isDead(theDetid, col, row) &&
0310 !theSiPixelGainCalibrationService_->isNoisy(theDetid, col, row)) {
0311
0312
0313
0314
0315
0316
0317
0318
0319 float DBgain = theSiPixelGainCalibrationService_->getGain(theDetid, col, row);
0320 float pedestal = theSiPixelGainCalibrationService_->getPedestal(theDetid, col, row);
0321 float DBpedestal = pedestal * DBgain;
0322
0323
0324
0325
0326
0327 float vcal = adc * DBgain - DBpedestal;
0328
0329
0330
0331
0332
0333
0334
0335
0336
0337
0338
0339
0340
0341
0342 if (theLayer == 1) {
0343 electrons = int(vcal * theConversionFactor_L1 + theOffset_L1);
0344 } else {
0345 electrons = int(vcal * theConversionFactor + theOffset);
0346 }
0347 }
0348 } else {
0349
0350 const float gain = theElectronPerADCGain;
0351 const float pedestal = 0.;
0352 electrons = int(adc * gain + pedestal);
0353 }
0354
0355 return electrons;
0356 }
0357
0358
0359
0360
0361 SiPixelCluster PixelThresholdClusterizer::make_cluster(const SiPixelCluster::PixelPos& pix,
0362 edmNew::DetSetVector<SiPixelCluster>::FastFiller& output) {
0363
0364 uint16_t seed_adc;
0365 std::stack<SiPixelCluster::PixelPos, std::vector<SiPixelCluster::PixelPos> > dead_pixel_stack;
0366
0367
0368
0369
0370
0371
0372
0373
0374
0375
0376
0377
0378
0379
0380
0381
0382
0383
0384
0385
0386
0387 seed_adc = std::min(theBuffer(pix.row(), pix.col()), int(std::numeric_limits<uint16_t>::max()));
0388 theBuffer.set_adc(pix, 1);
0389
0390
0391 AccretionCluster acluster, cldata;
0392 acluster.add(pix, seed_adc);
0393 if (!theFakePixels[pix.row() * theNumOfCols + pix.col()]) {
0394 cldata.add(pix, seed_adc);
0395 }
0396
0397
0398 bool dead_flag = false;
0399 while (!acluster.empty()) {
0400
0401 auto curInd = acluster.top();
0402 acluster.pop();
0403 for (auto c = std::max(0, int(acluster.y[curInd]) - 1);
0404 c < std::min(int(acluster.y[curInd]) + 2, theBuffer.columns());
0405 ++c) {
0406 for (auto r = std::max(0, int(acluster.x[curInd]) - 1);
0407 r < std::min(int(acluster.x[curInd]) + 2, theBuffer.rows());
0408 ++r) {
0409 if (theBuffer(r, c) >= thePixelThreshold) {
0410 SiPixelCluster::PixelPos newpix(r, c);
0411 auto const newpix_adc = std::min(theBuffer(r, c), int(std::numeric_limits<uint16_t>::max()));
0412 if (!acluster.add(newpix, newpix_adc))
0413 goto endClus;
0414
0415 if (!theFakePixels[r * theNumOfCols + c]) {
0416 cldata.add(newpix, newpix_adc);
0417 }
0418 theBuffer.set_adc(newpix, 1);
0419 }
0420
0421
0422
0423
0424
0425
0426
0427
0428
0429
0430
0431
0432
0433
0434
0435
0436
0437
0438
0439
0440
0441
0442
0443
0444 }
0445 }
0446
0447 }
0448 endClus:
0449 SiPixelCluster cluster(cldata.isize, cldata.adc, cldata.x, cldata.y, cldata.xmin, cldata.ymin);
0450
0451
0452 if (dead_flag && doSplitClusters) {
0453
0454 auto clusterThreshold = theClusterThreshold;
0455 if (theLayer == 1)
0456 clusterThreshold = theClusterThreshold_L1;
0457
0458
0459 SiPixelCluster first_cluster = cluster;
0460 bool have_second_cluster = false;
0461 while (!dead_pixel_stack.empty()) {
0462
0463 SiPixelCluster::PixelPos deadpix = dead_pixel_stack.top();
0464 dead_pixel_stack.pop();
0465 theBuffer.set_adc(deadpix, 1);
0466
0467
0468 SiPixelCluster second_cluster = make_cluster(deadpix, output);
0469
0470
0471 if (second_cluster.charge() >= clusterThreshold && first_cluster.charge() >= clusterThreshold) {
0472 output.push_back(second_cluster);
0473 have_second_cluster = true;
0474 }
0475
0476
0477
0478 const std::vector<SiPixelCluster::Pixel>& branch_pixels = second_cluster.pixels();
0479 for (unsigned int i = 0; i < branch_pixels.size(); i++) {
0480 auto const temp_x = branch_pixels[i].x;
0481 auto const temp_y = branch_pixels[i].y;
0482 auto const temp_adc = branch_pixels[i].adc;
0483 SiPixelCluster::PixelPos newpix(temp_x, temp_y);
0484 cluster.add(newpix, temp_adc);
0485 }
0486 }
0487
0488
0489 if (first_cluster.charge() >= clusterThreshold && have_second_cluster) {
0490 output.push_back(first_cluster);
0491 std::push_heap(output.begin(), output.end(), [](SiPixelCluster const& cl1, SiPixelCluster const& cl2) {
0492 return cl1.minPixelRow() < cl2.minPixelRow();
0493 });
0494 }
0495 }
0496
0497 return cluster;
0498 }