File indexing completed on 2023-03-17 11:19:34
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
0125
0126
0127
0128
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
0140 if (begin == end) {
0141 edm::LogError("PixelThresholdClusterizer") << "@SUB=PixelThresholdClusterizer::clusterizeDetUnitT()"
0142 << " No digis to clusterize";
0143 }
0144
0145
0146 if (!setup(pixDet))
0147 return;
0148
0149 theDetid = input.detId();
0150
0151
0152 auto clusterThreshold = theClusterThreshold;
0153 theLayer = (DetId(theDetid).subdetId() == 1) ? tTopo->pxbLayer(theDetid) : 0;
0154 if (theLayer == 1)
0155 clusterThreshold = theClusterThreshold_L1;
0156
0157
0158
0159 if (end > begin)
0160 copy_to_buffer(begin, end);
0161
0162 assert(output.empty());
0163
0164 for (unsigned int i = 0; i < theSeeds.size(); i++) {
0165
0166
0167 if (theBuffer(theSeeds[i]) >= theSeedThreshold) {
0168
0169 SiPixelCluster&& cluster = make_cluster(theSeeds[i], output);
0170
0171
0172
0173 if (cluster.charge() >= clusterThreshold) {
0174
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
0183 std::sort_heap(output.begin(), output.end(), [](SiPixelCluster const& cl1, SiPixelCluster const& cl2) {
0184 return cl1.minPixelRow() < cl2.minPixelRow();
0185 });
0186
0187
0188 theSeeds.clear();
0189
0190
0191 clear_buffer(begin, end);
0192
0193 theFakePixels.clear();
0194
0195 thePixelOccurrence.clear();
0196 }
0197
0198
0199
0200
0201
0202
0203
0204
0205
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);
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);
0219 }
0220 }
0221 }
0222
0223
0224
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
0232 }
0233 #endif
0234
0235
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];
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;
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
0282 int adc = (di->flag() != 0) ? di->adc() * 10 : electron[i];
0283 i++;
0284
0285 #ifdef PIXELREGRESSION
0286 int adcOld = calibrate(di->adc(), col, row);
0287
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;
0297
0298
0299
0300
0301 thePixelOccurrence[theBuffer.index(row, col)]++;
0302 uint8_t occurrence =
0303 (!dropDuplicates) ? 1 : thePixelOccurrence[theBuffer.index(row, col)];
0304
0305 switch (occurrence) {
0306
0307 case 1:
0308 if (adc >= thePixelThreshold) {
0309 theBuffer.set_adc(row, col, adc);
0310
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
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
0325 }
0326 }
0327 assert(i == (end - begin));
0328 }
0329
0330 void PixelThresholdClusterizer::copy_to_buffer(ClusterIterator begin, ClusterIterator end) {
0331
0332 for (ClusterIterator ci = begin; ci != end; ++ci) {
0333
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
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
0380
0381 if (!theSiPixelGainCalibrationService_->isDead(theDetid, col, row) &&
0382 !theSiPixelGainCalibrationService_->isNoisy(theDetid, col, row)) {
0383
0384
0385
0386
0387
0388
0389
0390
0391 float DBgain = theSiPixelGainCalibrationService_->getGain(theDetid, col, row);
0392 float pedestal = theSiPixelGainCalibrationService_->getPedestal(theDetid, col, row);
0393 float DBpedestal = pedestal * DBgain;
0394
0395
0396
0397
0398
0399 float vcal = adc * DBgain - DBpedestal;
0400
0401
0402
0403
0404
0405
0406
0407
0408
0409
0410
0411
0412
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 {
0421
0422 const float gain = theElectronPerADCGain;
0423 const float pedestal = 0.;
0424 electrons = int(adc * gain + pedestal);
0425 }
0426
0427 return electrons;
0428 }
0429
0430
0431
0432
0433 SiPixelCluster PixelThresholdClusterizer::make_cluster(const SiPixelCluster::PixelPos& pix,
0434 edmNew::DetSetVector<SiPixelCluster>::FastFiller& output) {
0435
0436 uint16_t seed_adc;
0437 std::stack<SiPixelCluster::PixelPos, std::vector<SiPixelCluster::PixelPos> > dead_pixel_stack;
0438
0439
0440
0441
0442
0443
0444
0445
0446
0447
0448
0449
0450
0451
0452
0453
0454
0455
0456
0457
0458
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
0468 bool dead_flag = false;
0469 while (!acluster.empty()) {
0470
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
0485 if (!theFakePixels[r * theNumOfCols + c]) {
0486 cldata.add(newpix, newpix_adc);
0487 }
0488 theBuffer.set_adc(newpix, 1);
0489 }
0490
0491
0492
0493
0494
0495
0496
0497
0498
0499
0500
0501
0502
0503
0504
0505
0506
0507
0508
0509
0510
0511
0512
0513
0514 }
0515 }
0516
0517 }
0518 endClus:
0519 SiPixelCluster cluster(cldata.isize, cldata.adc, cldata.x, cldata.y, cldata.xmin, cldata.ymin);
0520
0521
0522 if (dead_flag && doSplitClusters) {
0523
0524 auto clusterThreshold = theClusterThreshold;
0525 if (theLayer == 1)
0526 clusterThreshold = theClusterThreshold_L1;
0527
0528
0529 SiPixelCluster first_cluster = cluster;
0530 bool have_second_cluster = false;
0531 while (!dead_pixel_stack.empty()) {
0532
0533 SiPixelCluster::PixelPos deadpix = dead_pixel_stack.top();
0534 dead_pixel_stack.pop();
0535 theBuffer.set_adc(deadpix, 1);
0536
0537
0538 SiPixelCluster second_cluster = make_cluster(deadpix, output);
0539
0540
0541 if (second_cluster.charge() >= clusterThreshold && first_cluster.charge() >= clusterThreshold) {
0542 output.push_back(second_cluster);
0543 have_second_cluster = true;
0544 }
0545
0546
0547
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
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 }