File indexing completed on 2022-05-31 22:26:01
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 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
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
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
0092
0093
0094 bool PixelThresholdClusterizer::setup(const PixelGeomDetUnit* pixDet) {
0095
0096 const PixelTopology& topol = pixDet->specificTopology();
0097
0098
0099 int nrows = topol.nrows();
0100 int ncols = topol.ncolumns();
0101
0102 theNumOfRows = nrows;
0103 theNumOfCols = ncols;
0104
0105 if (nrows > theBuffer.rows() || ncols > theBuffer.columns()) {
0106 if (nrows != theNumOfRows || ncols != theNumOfCols)
0107 edm::LogWarning("setup()") << "pixel buffer redefined to" << nrows << " * " << ncols;
0108
0109
0110
0111 theBuffer.setSize(nrows, ncols);
0112 }
0113
0114 theFakePixels.resize(nrows * ncols, false);
0115
0116 return true;
0117 }
0118
0119
0120
0121
0122
0123
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
0135 if (begin == end) {
0136 edm::LogError("PixelThresholdClusterizer") << "@SUB=PixelThresholdClusterizer::clusterizeDetUnitT()"
0137 << " No digis to clusterize";
0138 }
0139
0140
0141 if (!setup(pixDet))
0142 return;
0143
0144 theDetid = input.detId();
0145
0146
0147 auto clusterThreshold = theClusterThreshold;
0148 theLayer = (DetId(theDetid).subdetId() == 1) ? tTopo->pxbLayer(theDetid) : 0;
0149 if (theLayer == 1)
0150 clusterThreshold = theClusterThreshold_L1;
0151
0152
0153
0154 if (end > begin)
0155 copy_to_buffer(begin, end);
0156
0157 assert(output.empty());
0158
0159 for (unsigned int i = 0; i < theSeeds.size(); i++) {
0160
0161
0162 if (theBuffer(theSeeds[i]) >= theSeedThreshold) {
0163
0164 SiPixelCluster&& cluster = make_cluster(theSeeds[i], output);
0165
0166
0167
0168 if (cluster.charge() >= clusterThreshold) {
0169
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
0178 std::sort_heap(output.begin(), output.end(), [](SiPixelCluster const& cl1, SiPixelCluster const& cl2) {
0179 return cl1.minPixelRow() < cl2.minPixelRow();
0180 });
0181
0182
0183 theSeeds.clear();
0184
0185
0186 clear_buffer(begin, end);
0187
0188 theFakePixels.clear();
0189 }
0190
0191
0192
0193
0194
0195
0196
0197
0198
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);
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);
0212 }
0213 }
0214 }
0215
0216
0217
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
0225 }
0226 #endif
0227
0228
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];
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;
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
0275 int adc = (di->flag() != 0) ? di->adc() * 10 : electron[i];
0276 i++;
0277
0278 #ifdef PIXELREGRESSION
0279 int adcOld = calibrate(di->adc(), col, row);
0280
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;
0290
0291
0292
0293
0294 if (adc >= thePixelThreshold) {
0295 theBuffer.set_adc(row, col, adc);
0296
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
0308 for (ClusterIterator ci = begin; ci != end; ++ci) {
0309
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
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
0356
0357 if (!theSiPixelGainCalibrationService_->isDead(theDetid, col, row) &&
0358 !theSiPixelGainCalibrationService_->isNoisy(theDetid, col, row)) {
0359
0360
0361
0362
0363
0364
0365
0366
0367 float DBgain = theSiPixelGainCalibrationService_->getGain(theDetid, col, row);
0368 float pedestal = theSiPixelGainCalibrationService_->getPedestal(theDetid, col, row);
0369 float DBpedestal = pedestal * DBgain;
0370
0371
0372
0373
0374
0375 float vcal = adc * DBgain - DBpedestal;
0376
0377
0378
0379
0380
0381
0382
0383
0384
0385
0386
0387
0388
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 {
0397
0398 const float gain = theElectronPerADCGain;
0399 const float pedestal = 0.;
0400 electrons = int(adc * gain + pedestal);
0401 }
0402
0403 return electrons;
0404 }
0405
0406
0407
0408
0409 SiPixelCluster PixelThresholdClusterizer::make_cluster(const SiPixelCluster::PixelPos& pix,
0410 edmNew::DetSetVector<SiPixelCluster>::FastFiller& output) {
0411
0412 uint16_t seed_adc;
0413 std::stack<SiPixelCluster::PixelPos, std::vector<SiPixelCluster::PixelPos> > dead_pixel_stack;
0414
0415
0416
0417
0418
0419
0420
0421
0422
0423
0424
0425
0426
0427
0428
0429
0430
0431
0432
0433
0434
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
0444 bool dead_flag = false;
0445 while (!acluster.empty()) {
0446
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
0461 if (!theFakePixels[r * theNumOfCols + c]) {
0462 cldata.add(newpix, newpix_adc);
0463 }
0464 theBuffer.set_adc(newpix, 1);
0465 }
0466
0467
0468
0469
0470
0471
0472
0473
0474
0475
0476
0477
0478
0479
0480
0481
0482
0483
0484
0485
0486
0487
0488
0489
0490 }
0491 }
0492
0493 }
0494 endClus:
0495 SiPixelCluster cluster(cldata.isize, cldata.adc, cldata.x, cldata.y, cldata.xmin, cldata.ymin);
0496
0497
0498 if (dead_flag && doSplitClusters) {
0499
0500 auto clusterThreshold = theClusterThreshold;
0501 if (theLayer == 1)
0502 clusterThreshold = theClusterThreshold_L1;
0503
0504
0505 SiPixelCluster first_cluster = cluster;
0506 bool have_second_cluster = false;
0507 while (!dead_pixel_stack.empty()) {
0508
0509 SiPixelCluster::PixelPos deadpix = dead_pixel_stack.top();
0510 dead_pixel_stack.pop();
0511 theBuffer.set_adc(deadpix, 1);
0512
0513
0514 SiPixelCluster second_cluster = make_cluster(deadpix, output);
0515
0516
0517 if (second_cluster.charge() >= clusterThreshold && first_cluster.charge() >= clusterThreshold) {
0518 output.push_back(second_cluster);
0519 have_second_cluster = true;
0520 }
0521
0522
0523
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
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 }