Skip to content

Commit

Permalink
Backport of changes in PR34606 merged into 12_1_X
Browse files Browse the repository at this point in the history
  • Loading branch information
veszpv committed Aug 10, 2021
1 parent 11f2a69 commit db8c328
Show file tree
Hide file tree
Showing 9 changed files with 262 additions and 4 deletions.
@@ -0,0 +1,4 @@
import FWCore.ParameterSet.Config as cms

# This modifier is to run the SiPixelDigiMorphing
siPixelDigiMorphing = cms.Modifier()
4 changes: 4 additions & 0 deletions RecoLocalTracker/Configuration/python/RecoLocalTracker_cff.py
Expand Up @@ -8,13 +8,17 @@
from RecoLocalTracker.SiStripZeroSuppression.SiStripZeroSuppression_cfi import *
from RecoLocalTracker.SiStripClusterizer.SiStripClusterizer_cfi import *
from RecoLocalTracker.SiPixelClusterizer.siPixelClustersPreSplitting_cff import *
from RecoLocalTracker.SiPixelDigiReProducers.siPixelDigisMorphed_cfi import *
from RecoLocalTracker.SiPixelRecHits.SiPixelRecHits_cfi import *
from RecoLocalTracker.SubCollectionProducers.clustersummaryproducer_cfi import *

pixeltrackerlocalrecoTask = cms.Task(
siPixelClustersPreSplittingTask,
siPixelRecHitsPreSplittingTask)

from Configuration.ProcessModifiers.siPixelDigiMorphing_cff import *
siPixelDigiMorphing.toModify(pixeltrackerlocalrecoTask, func=lambda t: t.add(siPixelDigisMorphed))

striptrackerlocalrecoTask = cms.Task(
siStripZeroSuppression,
siStripClusters,
Expand Down
Expand Up @@ -62,6 +62,7 @@ PixelThresholdClusterizer::PixelThresholdClusterizer(edm::ParameterSet const& co
doMissCalibrate(conf.getParameter<bool>("MissCalibrate")),
doSplitClusters(conf.getParameter<bool>("SplitClusters")) {
theBuffer.setSize(theNumOfRows, theNumOfCols);
theFakePixels.clear();
}
/////////////////////////////////////////////////////////////////////////////
PixelThresholdClusterizer::~PixelThresholdClusterizer() {}
Expand Down Expand Up @@ -110,6 +111,8 @@ bool PixelThresholdClusterizer::setup(const PixelGeomDetUnit* pixDet) {
theBuffer.setSize(nrows, ncols); // Modify
}

theFakePixels.resize(nrows * ncols, false);

return true;
}
//----------------------------------------------------------------------------
Expand Down Expand Up @@ -179,6 +182,8 @@ void PixelThresholdClusterizer::clusterizeDetUnitT(const T& input,

// Need to clean unused pixels from the buffer array.
clear_buffer(begin, end);

theFakePixels.clear();
}

//----------------------------------------------------------------------------
Expand Down Expand Up @@ -218,7 +223,7 @@ void PixelThresholdClusterizer::copy_to_buffer(DigiIterator begin, DigiIterator
}
#endif
int electron[end - begin]; // pixel charge in electrons
memset(electron, 0, sizeof(electron));
memset(electron, 0, (end - begin) * sizeof(int));

if (doPhase2Calibration) {
int i = 0;
Expand Down Expand Up @@ -257,7 +262,9 @@ void PixelThresholdClusterizer::copy_to_buffer(DigiIterator begin, DigiIterator
for (DigiIterator di = begin; di != end; ++di) {
int row = di->row();
int col = di->column();
int adc = electron[i++]; // this is in electrons
// VV: do not calibrate a fake pixel, it already has a unit of 10e-:
int adc = (di->flag() != 0) ? di->adc() * 10 : electron[i]; // this is in electrons
i++;

#ifdef PIXELREGRESSION
int adcOld = calibrate(di->adc(), col, row);
Expand All @@ -277,6 +284,9 @@ void PixelThresholdClusterizer::copy_to_buffer(DigiIterator begin, DigiIterator

if (adc >= thePixelThreshold) {
theBuffer.set_adc(row, col, adc);
// VV: add pixel to the fake list. Only when running on digi collection
if (di->flag() != 0)
theFakePixels[row * theNumOfCols + col] = true;
if (adc >= theSeedThreshold)
theSeeds.push_back(SiPixelCluster::PixelPos(row, col));
}
Expand Down Expand Up @@ -414,8 +424,9 @@ SiPixelCluster PixelThresholdClusterizer::make_cluster(const SiPixelCluster::Pix
theBuffer.set_adc(pix, 1);
// }

AccretionCluster acluster;
AccretionCluster acluster, cldata;
acluster.add(pix, seed_adc);
cldata.add(pix, seed_adc);

//Here we search all pixels adjacent to all pixels in the cluster.
bool dead_flag = false;
Expand All @@ -433,6 +444,10 @@ SiPixelCluster PixelThresholdClusterizer::make_cluster(const SiPixelCluster::Pix
SiPixelCluster::PixelPos newpix(r, c);
if (!acluster.add(newpix, theBuffer(r, c)))
goto endClus;
// VV: no fake pixels in cluster, leads to non-contiguous clusters
if (!theFakePixels[r * theNumOfCols + c]) {
cldata.add(newpix, theBuffer(r, c));
}
theBuffer.set_adc(newpix, 1);
}

Expand Down Expand Up @@ -464,7 +479,7 @@ SiPixelCluster PixelThresholdClusterizer::make_cluster(const SiPixelCluster::Pix

} // while accretion
endClus:
SiPixelCluster cluster(acluster.isize, acluster.adc, acluster.x, acluster.y, acluster.xmin, acluster.ymin);
SiPixelCluster cluster(cldata.isize, cldata.adc, cldata.x, cldata.y, cldata.xmin, cldata.ymin);
//Here we split the cluster, if the flag to do so is set and we have found a dead or noisy pixel.

if (dead_flag && doSplitClusters) {
Expand Down
Expand Up @@ -90,6 +90,8 @@ class dso_hidden PixelThresholdClusterizer final : public PixelClusterizerBase {
std::vector<SiPixelCluster::PixelPos> theSeeds; // cached seed pixels
std::vector<SiPixelCluster> theClusters; // resulting clusters

std::vector<bool> theFakePixels; // fake pixels introduced to guide clustering

//! Clustering-related quantities:
float thePixelThresholdInNoiseUnits; // Pixel threshold in units of noise
float theSeedThresholdInNoiseUnits; // Pixel cluster seed in units of noise
Expand Down
Expand Up @@ -23,3 +23,11 @@
)
)
)

from Configuration.ProcessModifiers.siPixelDigiMorphing_cff import siPixelDigiMorphing
siPixelDigiMorphing.toModify(
siPixelClustersPreSplitting,
cpu = dict(
src = 'siPixelDigisMorphed'
)
)
@@ -0,0 +1,6 @@

/*!

\page RecoLocalTracker_SiPixelDigiReProducers Package RecoLocalTracker/SiPixelDigiReProducers
*/

7 changes: 7 additions & 0 deletions RecoLocalTracker/SiPixelDigiReProducers/plugins/BuildFile.xml
@@ -0,0 +1,7 @@
<use name="DataFormats/Common"/>
<use name="FWCore/ParameterSet"/>
<use name="DataFormats/SiPixelDigi"/>
<use name="CalibTracker/SiPixelESProducers"/>
<library file="*.cc" name="RecoLocalTrackerSiPixelDigiReProducersPlugins">
<flags EDM_PLUGIN="1"/>
</library>
207 changes: 207 additions & 0 deletions RecoLocalTracker/SiPixelDigiReProducers/plugins/SiPixelDigiMorphing.cc
@@ -0,0 +1,207 @@
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/stream/EDProducer.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/PluginManager/interface/ModuleDef.h"
#include "DataFormats/Common/interface/DetSetVector.h"
#include "DataFormats/Common/interface/DetSetVectorNew.h"
#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/SiPixelDigi/interface/PixelDigi.h"

#include <bitset>
#include <iterator>

class SiPixelDigiMorphing : public edm::stream::EDProducer<> {
public:
explicit SiPixelDigiMorphing(const edm::ParameterSet& conf);
~SiPixelDigiMorphing() override = default;

static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
void produce(edm::Event& e, const edm::EventSetup& c) override;

private:
edm::EDGetTokenT<edm::DetSetVector<PixelDigi>> tPixelDigi_;
edm::EDPutTokenT<edm::DetSetVector<PixelDigi>> tPutPixelDigi_;

const int32_t nrows_;
const int32_t ncols_;
const int32_t nrocs_; // in Phase 1, this is ROCs, but could be any subset of a pixel row
const int32_t iters_;
const uint32_t fakeAdc_;

int32_t ncols_r_; // number of columns per ROC
int32_t ksize_; // kernel size

std::vector<uint64_t> kernel1_;
std::vector<uint64_t> kernel2_;
uint64_t mask_;

enum MorphOption { kDilate, kErode };

void morph(uint64_t* const imap, uint64_t* omap, uint64_t* const kernel, MorphOption op) const;
};

SiPixelDigiMorphing::SiPixelDigiMorphing(edm::ParameterSet const& conf)
: tPixelDigi_(consumes(conf.getParameter<edm::InputTag>("src"))),
tPutPixelDigi_(produces<edm::DetSetVector<PixelDigi>>()),
nrows_(conf.getParameter<int32_t>("nrows")),
ncols_(conf.getParameter<int32_t>("ncols")),
nrocs_(conf.getParameter<int32_t>("nrocs")),
iters_(conf.getParameter<int32_t>("iters")),
fakeAdc_(conf.getParameter<uint32_t>("fakeAdc")) {
if (ncols_ % nrocs_ == 0) {
ncols_r_ = ncols_ / nrocs_;
} else {
throw cms::Exception("Configuration") << "[SiPixelDigiMorphing]:"
<< " number of columns not divisible with"
<< " number of ROCs\n";
}

if (ncols_r_ + 2 * iters_ <= int(sizeof(uint64_t) * 8)) {
ksize_ = 2 * iters_ + 1;
} else {
throw cms::Exception("Configuration") << "[SiPixelDigiMorphing]:"
<< " too many columns per ROC"
<< " or too many iterations set\n"
<< " Ncol/Nrocs+2*iters should not be"
<< " more than " << sizeof(uint64_t) * 8 << "\n";
}

std::vector<int32_t> k1(conf.getParameter<std::vector<int32_t>>("kernel1"));
std::vector<int32_t> k2(conf.getParameter<std::vector<int32_t>>("kernel2"));

kernel1_.resize(ksize_, 0);
kernel2_.resize(ksize_, 0);
mask_ = 0;
int w = (ncols_r_ + 2 * iters_) / ksize_ + ((ncols_r_ + 2 * iters_) % ksize_ != 0);
for (int j = 0; j < w; j++) {
for (int ii = 0; ii < ksize_; ii++) {
kernel1_[ii] <<= ksize_;
kernel1_[ii] |= k1[ii];
kernel2_[ii] <<= ksize_;
kernel2_[ii] |= k2[ii];
}
mask_ <<= ksize_;
mask_ |= 1;
}
mask_ <<= iters_;
}

void SiPixelDigiMorphing::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
edm::ParameterSetDescription desc;

desc.add<edm::InputTag>("src", edm::InputTag("siPixelDigis"));
desc.add<int32_t>("nrows", 160);
desc.add<int32_t>("ncols", 416);
desc.add<int32_t>("nrocs", 8);
desc.add<int32_t>("iters", 1);
desc.add<std::vector<int32_t>>("kernel1", {7, 7, 7});
desc.add<std::vector<int32_t>>("kernel2", {2, 7, 2});
desc.add<uint32_t>("fakeAdc", 100);

descriptions.addWithDefaultLabel(desc);
}

void SiPixelDigiMorphing::produce(edm::Event& e, const edm::EventSetup& es) {
auto const& inputDigi = e.get(tPixelDigi_);

auto outputDigis = std::make_unique<edm::DetSetVector<PixelDigi>>();

const int rocSize = nrows_ + 2 * iters_;
const int arrSize = nrocs_ * rocSize;

uint64_t imap[arrSize];
uint64_t map1[arrSize];
uint64_t map2[arrSize];

for (auto const& ds : inputDigi) {
auto rawId = ds.detId();
edm::DetSet<PixelDigi>* detDigis = nullptr;
detDigis = &(outputDigis->find_or_insert(rawId));

memset(imap, 0, arrSize * sizeof(uint64_t));
for (auto const& di : ds) {
int r = int(di.column()) / ncols_r_;
int c = int(di.column()) % ncols_r_;
imap[r * rocSize + di.row() + iters_] |= uint64_t(1) << (c + iters_);
if (r > 0 && c < iters_) {
imap[(r - 1) * rocSize + di.row() + iters_] |= uint64_t(1) << (c + ncols_r_ + iters_);
} else if (++r < nrocs_ && c >= ncols_r_ - iters_) {
imap[r * rocSize + di.row() + iters_] |= uint64_t(1) << (c - ncols_r_ + iters_);
}
(*detDigis).data.emplace_back(di.row(), di.column(), di.adc(), 0);
}

std::memcpy(map1, imap, arrSize * sizeof(uint64_t));
memset(map2, 0, arrSize * sizeof(uint64_t));

morph(map1, map2, kernel1_.data(), kDilate);
morph(map2, map1, kernel2_.data(), kErode);

uint64_t* i = imap + iters_;
uint64_t* o = map1 + iters_;
for (int roc = 0; roc < nrocs_; roc++, i += 2 * iters_, o += 2 * iters_) {
for (int row = 0; row < nrows_; row++, i++, o++) {
if (*o == 0)
continue;
*o >>= iters_;
*i >>= iters_;
for (int col = 0; col < ncols_r_; col++, (*i) >>= 1, (*o) >>= 1) {
if (*o == 0)
break;
if (((*i) & uint64_t(1)) == 1)
continue;
if (((*o) & uint64_t(1)) == 1) {
(*detDigis).data.emplace_back(row, roc * ncols_r_ + col, fakeAdc_, 1);
}
}
}
}
}
e.put(tPutPixelDigi_, std::move(outputDigis));
}

void SiPixelDigiMorphing::morph(uint64_t* const imap, uint64_t* omap, uint64_t* const kernel, MorphOption op) const {
uint64_t* i[ksize_]; // i(nput)
uint64_t* o = omap + iters_; // o(output)
unsigned char valid = 0;
unsigned char const validMask = (1 << ksize_) - 1;
uint64_t m[ksize_]; // m(ask)

for (int ii = 0; ii < ksize_; ii++) {
i[ii] = imap + ii;
valid = (valid << 1) | (*i[ii] != 0);
m[ii] = mask_ << ii;
}

for (int roc = 0; roc < nrocs_; roc++, o += 2 * iters_) {
for (int row = 0; row < nrows_; row++, o++) {
if ((valid & validMask) != 0) {
for (int jj = 0; jj < ksize_; jj++) {
for (int ii = 0; ii < ksize_; ii++) {
uint64_t v = (*i[ii]) & (kernel[ii] << jj); // v(ector)
if (op == kErode)
v ^= (kernel[ii] << jj);
uint64_t vv = v; // vv(vector bit - shifted and contracted)
for (int b = 1; b < ksize_; b++)
vv |= (v >> b);
*o |= ((vv << iters_) & m[jj]);
}
if (op == kErode)
*o ^= m[jj];
}
}
for (int ii = 0; ii < ksize_; ii++)
i[ii]++;
valid = (valid << 1) | (*i[ksize_ - 1] != 0);
}
for (int ii = 0; ii < ksize_; ii++) {
i[ii] += 2 * iters_;
valid = (valid << 1) | (*i[ii] != 0);
}
}
}

DEFINE_FWK_MODULE(SiPixelDigiMorphing);
@@ -0,0 +1,5 @@
import FWCore.ParameterSet.Config as cms

from RecoLocalTracker.SiPixelDigiReProducers.siPixelDigiMorphing_cfi import siPixelDigiMorphing as _siPixelDigiMorphing

siPixelDigisMorphed = _siPixelDigiMorphing.clone()

0 comments on commit db8c328

Please sign in to comment.