File indexing completed on 2023-03-17 11:19:52
0001 #include "FWCore/Framework/interface/stream/EDProducer.h"
0002
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005 #include "FWCore/Utilities/interface/InputTag.h"
0006 #include "DataFormats/Common/interface/Handle.h"
0007 #include "FWCore/Framework/interface/ESHandle.h"
0008 #include "DataFormats/SiPixelCluster/interface/SiPixelCluster.h"
0009 #include "DataFormats/Common/interface/DetSetVectorNew.h"
0010
0011 #include "RecoLocalTracker/ClusterParameterEstimator/interface/PixelClusterParameterEstimator.h"
0012 #include "RecoLocalTracker/Records/interface/TkPixelCPERecord.h"
0013
0014 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0015 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0016 #include "DataFormats/GeometryVector/interface/VectorUtil.h"
0017
0018 #include "DataFormats/VertexReco/interface/Vertex.h"
0019 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0020 #include "DataFormats/JetReco/interface/Jet.h"
0021
0022 #include <algorithm>
0023 #include <vector>
0024 #include <utility>
0025
0026 class JetCoreClusterSplitter : public edm::stream::EDProducer<> {
0027 public:
0028 JetCoreClusterSplitter(const edm::ParameterSet& iConfig);
0029 ~JetCoreClusterSplitter() override;
0030 void produce(edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0031
0032 private:
0033 bool split(const SiPixelCluster& aCluster,
0034 edmNew::DetSetVector<SiPixelCluster>::FastFiller& filler,
0035 float expectedADC,
0036 int sizeY,
0037 int sizeX,
0038 float jetZOverRho);
0039 std::vector<SiPixelCluster> fittingSplit(const SiPixelCluster& aCluster,
0040 float expectedADC,
0041 int sizeY,
0042 int sizeX,
0043 float jetZOverRho,
0044 unsigned int nSplitted);
0045 std::pair<float, float> closestClusters(const std::vector<float>& distanceMap);
0046 std::multimap<float, int> secondDistDiffScore(const std::vector<std::vector<float>>& distanceMap);
0047 std::multimap<float, int> secondDistScore(const std::vector<std::vector<float>>& distanceMap);
0048 std::multimap<float, int> distScore(const std::vector<std::vector<float>>& distanceMap);
0049
0050 edm::ESGetToken<GlobalTrackingGeometry, GlobalTrackingGeometryRecord> const tTrackingGeom_;
0051 edm::ESGetToken<PixelClusterParameterEstimator, TkPixelCPERecord> const tCPE_;
0052
0053 bool verbose;
0054 double ptMin_;
0055 double deltaR_;
0056 double chargeFracMin_;
0057 edm::EDGetTokenT<edmNew::DetSetVector<SiPixelCluster>> pixelClusters_;
0058 edm::EDGetTokenT<reco::VertexCollection> vertices_;
0059 edm::EDGetTokenT<edm::View<reco::Candidate>> cores_;
0060 double forceXError_;
0061 double forceYError_;
0062 double fractionalWidth_;
0063 double chargePerUnit_;
0064 double centralMIPCharge_;
0065 };
0066
0067 JetCoreClusterSplitter::JetCoreClusterSplitter(const edm::ParameterSet& iConfig)
0068 : tTrackingGeom_(esConsumes()),
0069 tCPE_(esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("pixelCPE")))),
0070 verbose(iConfig.getParameter<bool>("verbose")),
0071 ptMin_(iConfig.getParameter<double>("ptMin")),
0072 deltaR_(iConfig.getParameter<double>("deltaRmax")),
0073 chargeFracMin_(iConfig.getParameter<double>("chargeFractionMin")),
0074 pixelClusters_(
0075 consumes<edmNew::DetSetVector<SiPixelCluster>>(iConfig.getParameter<edm::InputTag>("pixelClusters"))),
0076 vertices_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertices"))),
0077 cores_(consumes<edm::View<reco::Candidate>>(iConfig.getParameter<edm::InputTag>("cores"))),
0078 forceXError_(iConfig.getParameter<double>("forceXError")),
0079 forceYError_(iConfig.getParameter<double>("forceYError")),
0080 fractionalWidth_(iConfig.getParameter<double>("fractionalWidth")),
0081 chargePerUnit_(iConfig.getParameter<double>("chargePerUnit")),
0082 centralMIPCharge_(iConfig.getParameter<double>("centralMIPCharge"))
0083
0084 {
0085 produces<edmNew::DetSetVector<SiPixelCluster>>();
0086 }
0087
0088 JetCoreClusterSplitter::~JetCoreClusterSplitter() {}
0089
0090 bool SortPixels(const SiPixelCluster::Pixel& i, const SiPixelCluster::Pixel& j) { return (i.adc > j.adc); }
0091
0092 void JetCoreClusterSplitter::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0093 using namespace edm;
0094 const auto& geometry = &iSetup.getData(tTrackingGeom_);
0095
0096 Handle<edmNew::DetSetVector<SiPixelCluster>> inputPixelClusters;
0097 iEvent.getByToken(pixelClusters_, inputPixelClusters);
0098
0099 Handle<std::vector<reco::Vertex>> vertices;
0100 iEvent.getByToken(vertices_, vertices);
0101 const reco::Vertex& pv = (*vertices)[0];
0102
0103 Handle<edm::View<reco::Candidate>> cores;
0104 iEvent.getByToken(cores_, cores);
0105
0106 const PixelClusterParameterEstimator* pp = &iSetup.getData(tCPE_);
0107 auto output = std::make_unique<edmNew::DetSetVector<SiPixelCluster>>();
0108
0109 edmNew::DetSetVector<SiPixelCluster>::const_iterator detIt = inputPixelClusters->begin();
0110 for (; detIt != inputPixelClusters->end(); detIt++) {
0111 edmNew::DetSetVector<SiPixelCluster>::FastFiller filler(*output, detIt->id());
0112 const edmNew::DetSet<SiPixelCluster>& detset = *detIt;
0113 const GeomDet* det = geometry->idToDet(detset.id());
0114 for (auto cluster = detset.begin(); cluster != detset.end(); cluster++) {
0115 const SiPixelCluster& aCluster = *cluster;
0116 bool hasBeenSplit = false;
0117 bool shouldBeSplit = false;
0118 GlobalPoint cPos =
0119 det->surface().toGlobal(pp->localParametersV(aCluster, (*geometry->idToDetUnit(detIt->id())))[0].first);
0120 GlobalPoint ppv(pv.position().x(), pv.position().y(), pv.position().z());
0121 GlobalVector clusterDir = cPos - ppv;
0122 for (unsigned int ji = 0; ji < cores->size(); ji++) {
0123 if ((*cores)[ji].pt() > ptMin_) {
0124 const reco::Candidate& jet = (*cores)[ji];
0125 GlobalVector jetDir(jet.px(), jet.py(), jet.pz());
0126 if (Geom::deltaR(jetDir, clusterDir) < deltaR_) {
0127
0128
0129 bool isEndCap = (std::abs(cPos.z()) > 30.f);
0130 float jetZOverRho = jet.momentum().Z() / jet.momentum().Rho();
0131 if (isEndCap)
0132 jetZOverRho = jet.momentum().Rho() / jet.momentum().Z();
0133 float expSizeY = std::sqrt((1.3f * 1.3f) + (1.9f * 1.9f) * jetZOverRho * jetZOverRho);
0134 if (expSizeY < 1.f)
0135 expSizeY = 1.f;
0136 float expSizeX = 1.5f;
0137 if (isEndCap) {
0138 expSizeX = expSizeY;
0139 expSizeY = 1.5f;
0140 }
0141 float expCharge = std::sqrt(1.08f + jetZOverRho * jetZOverRho) * centralMIPCharge_;
0142
0143 if (aCluster.charge() > expCharge * chargeFracMin_ &&
0144 (aCluster.sizeX() > expSizeX + 1 || aCluster.sizeY() > expSizeY + 1)) {
0145 shouldBeSplit = true;
0146 if (verbose)
0147 std::cout << "Trying to split: charge and deltaR " << aCluster.charge() << " "
0148 << Geom::deltaR(jetDir, clusterDir) << " size x y " << aCluster.sizeX() << " "
0149 << aCluster.sizeY() << " exp. size (x,y) " << expSizeX << " " << expSizeY << " detid "
0150 << detIt->id() << std::endl;
0151 if (verbose)
0152 std::cout << "jetZOverRho=" << jetZOverRho << std::endl;
0153
0154 if (split(aCluster, filler, expCharge, expSizeY, expSizeX, jetZOverRho)) {
0155 hasBeenSplit = true;
0156 }
0157 }
0158 }
0159 }
0160 }
0161 if (!hasBeenSplit) {
0162 SiPixelCluster c = aCluster;
0163 if (shouldBeSplit) {
0164
0165
0166 c.setSplitClusterErrorX(c.sizeX() * (100.f / 3.f));
0167 c.setSplitClusterErrorY(c.sizeY() * (150.f / 3.f));
0168 }
0169 filler.push_back(c);
0170 std::push_heap(filler.begin(), filler.end(), [](SiPixelCluster const& cl1, SiPixelCluster const& cl2) {
0171 return cl1.minPixelRow() < cl2.minPixelRow();
0172 });
0173 }
0174 }
0175 std::sort_heap(filler.begin(), filler.end(), [](SiPixelCluster const& cl1, SiPixelCluster const& cl2) {
0176 return cl1.minPixelRow() < cl2.minPixelRow();
0177 });
0178
0179 }
0180 iEvent.put(std::move(output));
0181 }
0182
0183 bool JetCoreClusterSplitter::split(const SiPixelCluster& aCluster,
0184 edmNew::DetSetVector<SiPixelCluster>::FastFiller& filler,
0185 float expectedADC,
0186 int sizeY,
0187 int sizeX,
0188 float jetZOverRho) {
0189
0190
0191
0192 std::vector<SiPixelCluster> sp = fittingSplit(
0193 aCluster, expectedADC, sizeY, sizeX, jetZOverRho, std::floor(aCluster.charge() / expectedADC + 0.5f));
0194
0195
0196 for (unsigned int i = 0; i < sp.size(); i++) {
0197 filler.push_back(sp[i]);
0198 std::push_heap(filler.begin(), filler.end(), [](SiPixelCluster const& cl1, SiPixelCluster const& cl2) {
0199 return cl1.minPixelRow() < cl2.minPixelRow();
0200 });
0201 }
0202
0203 return (!sp.empty());
0204 }
0205
0206
0207 std::pair<float, float> JetCoreClusterSplitter::closestClusters(const std::vector<float>& distanceMap) {
0208 float minDist = std::numeric_limits<float>::max();
0209 float secondMinDist = std::numeric_limits<float>::max();
0210 for (unsigned int i = 0; i < distanceMap.size(); i++) {
0211 float dist = distanceMap[i];
0212 if (dist < minDist) {
0213 secondMinDist = minDist;
0214 minDist = dist;
0215 } else if (dist < secondMinDist) {
0216 secondMinDist = dist;
0217 }
0218 }
0219 return std::pair<float, float>(minDist, secondMinDist);
0220 }
0221
0222 std::multimap<float, int> JetCoreClusterSplitter::secondDistDiffScore(
0223 const std::vector<std::vector<float>>& distanceMap) {
0224 std::multimap<float, int> scores;
0225 for (unsigned int j = 0; j < distanceMap.size(); j++) {
0226 std::pair<float, float> d = closestClusters(distanceMap[j]);
0227 scores.insert(std::pair<float, int>(d.second - d.first, j));
0228 }
0229 return scores;
0230 }
0231
0232 std::multimap<float, int> JetCoreClusterSplitter::secondDistScore(const std::vector<std::vector<float>>& distanceMap) {
0233 std::multimap<float, int> scores;
0234 for (unsigned int j = 0; j < distanceMap.size(); j++) {
0235 std::pair<float, float> d = closestClusters(distanceMap[j]);
0236 scores.insert(std::pair<float, int>(-d.second, j));
0237 }
0238 return scores;
0239 }
0240
0241 std::multimap<float, int> JetCoreClusterSplitter::distScore(const std::vector<std::vector<float>>& distanceMap) {
0242 std::multimap<float, int> scores;
0243 for (unsigned int j = 0; j < distanceMap.size(); j++) {
0244 std::pair<float, float> d = closestClusters(distanceMap[j]);
0245 scores.insert(std::pair<float, int>(-d.first, j));
0246 }
0247 return scores;
0248 }
0249
0250 std::vector<SiPixelCluster> JetCoreClusterSplitter::fittingSplit(const SiPixelCluster& aCluster,
0251 float expectedADC,
0252 int sizeY,
0253 int sizeX,
0254 float jetZOverRho,
0255 unsigned int nSplitted) {
0256 std::vector<SiPixelCluster> output;
0257
0258 unsigned int meanExp = nSplitted;
0259 if (meanExp <= 1) {
0260 output.push_back(aCluster);
0261 return output;
0262 }
0263
0264 std::vector<float> clx(meanExp);
0265 std::vector<float> cly(meanExp);
0266 std::vector<float> cls(meanExp);
0267 std::vector<float> oldclx(meanExp);
0268 std::vector<float> oldcly(meanExp);
0269 std::vector<SiPixelCluster::Pixel> originalpixels = aCluster.pixels();
0270 std::vector<std::pair<int, SiPixelCluster::Pixel>> pixels;
0271 for (unsigned int j = 0; j < originalpixels.size(); j++) {
0272 int sub = originalpixels[j].adc / chargePerUnit_ * expectedADC / centralMIPCharge_;
0273 if (sub < 1)
0274 sub = 1;
0275 int perDiv = originalpixels[j].adc / sub;
0276 if (verbose)
0277 std::cout << "Splitting " << j << " in [ " << pixels.size() << " , " << pixels.size() + sub
0278 << " ], expected numb of clusters: " << meanExp << " original pixel (x,y) " << originalpixels[j].x
0279 << " " << originalpixels[j].y << " sub " << sub << std::endl;
0280 for (int k = 0; k < sub; k++) {
0281 if (k == sub - 1)
0282 perDiv = originalpixels[j].adc - perDiv * k;
0283 pixels.push_back(std::make_pair(j, SiPixelCluster::Pixel(originalpixels[j].x, originalpixels[j].y, perDiv)));
0284 }
0285 }
0286 std::vector<int> clusterForPixel(pixels.size());
0287
0288 for (unsigned int j = 0; j < meanExp; j++) {
0289 oldclx[j] = -999;
0290 oldcly[j] = -999;
0291 clx[j] = originalpixels[0].x + j;
0292 cly[j] = originalpixels[0].y + j;
0293 cls[j] = 0;
0294 }
0295 bool stop = false;
0296 int remainingSteps = 100;
0297 while (!stop && remainingSteps > 0) {
0298 remainingSteps--;
0299
0300 std::vector<std::vector<float>> distanceMapX(originalpixels.size(), std::vector<float>(meanExp));
0301 std::vector<std::vector<float>> distanceMapY(originalpixels.size(), std::vector<float>(meanExp));
0302 std::vector<std::vector<float>> distanceMap(originalpixels.size(), std::vector<float>(meanExp));
0303 for (unsigned int j = 0; j < originalpixels.size(); j++) {
0304 if (verbose)
0305 std::cout << "Original Pixel pos " << j << " " << pixels[j].second.x << " " << pixels[j].second.y << std::endl;
0306 for (unsigned int i = 0; i < meanExp; i++) {
0307 distanceMapX[j][i] = 1.f * originalpixels[j].x - clx[i];
0308 distanceMapY[j][i] = 1.f * originalpixels[j].y - cly[i];
0309 float dist = 0;
0310
0311 if (std::abs(distanceMapX[j][i]) > sizeX / 2.f) {
0312 dist +=
0313 (std::abs(distanceMapX[j][i]) - sizeX / 2.f + 1.f) * (std::abs(distanceMapX[j][i]) - sizeX / 2.f + 1.f);
0314 } else {
0315 dist += (2.f * distanceMapX[j][i] / sizeX) * (2.f * distanceMapX[j][i] / sizeX);
0316 }
0317
0318 if (std::abs(distanceMapY[j][i]) > sizeY / 2.f) {
0319 dist += 1.f * (std::abs(distanceMapY[j][i]) - sizeY / 2.f + 1.f) *
0320 (std::abs(distanceMapY[j][i]) - sizeY / 2.f + 1.f);
0321 } else {
0322 dist += 1.f * (2.f * distanceMapY[j][i] / sizeY) * (2.f * distanceMapY[j][i] / sizeY);
0323 }
0324 distanceMap[j][i] = sqrt(dist);
0325 if (verbose)
0326 std::cout << "Cluster " << i << " Original Pixel " << j << " distances: " << distanceMapX[j][i] << " "
0327 << distanceMapY[j][i] << " " << distanceMap[j][i] << std::endl;
0328 }
0329 }
0330
0331
0332
0333 std::multimap<float, int> scores;
0334
0335
0336 scores = secondDistScore(distanceMap);
0337
0338
0339
0340 std::vector<float> weightOfPixel(pixels.size());
0341 for (std::multimap<float, int>::iterator it = scores.begin(); it != scores.end(); it++) {
0342 int pixel_index = it->second;
0343 if (verbose)
0344 std::cout << "Original Pixel " << pixel_index << " with score " << it->first << std::endl;
0345
0346 int subpixel_counter = 0;
0347 for (auto subpixel = pixels.begin(); subpixel != pixels.end(); ++subpixel, ++subpixel_counter) {
0348 if (subpixel->first > pixel_index) {
0349 break;
0350 } else if (subpixel->first != pixel_index) {
0351 continue;
0352 } else {
0353 float maxEst = 0;
0354 int cl = -1;
0355 for (unsigned int subcluster_index = 0; subcluster_index < meanExp; subcluster_index++) {
0356 float nsig = (cls[subcluster_index] - expectedADC) /
0357 (expectedADC * fractionalWidth_);
0358 float clQest = 1.f / (1.f + std::exp(nsig)) + 1e-6f;
0359 float clDest = 1.f / (distanceMap[pixel_index][subcluster_index] + 0.05f);
0360
0361 if (verbose)
0362 std::cout << " Q: " << clQest << " D: " << clDest << " " << distanceMap[pixel_index][subcluster_index]
0363 << std::endl;
0364 float est = clQest * clDest;
0365 if (est > maxEst) {
0366 cl = subcluster_index;
0367 maxEst = est;
0368 }
0369 }
0370 cls[cl] += subpixel->second.adc;
0371 clusterForPixel[subpixel_counter] = cl;
0372 weightOfPixel[subpixel_counter] = maxEst;
0373 if (verbose)
0374 std::cout << "Pixel weight j cl " << weightOfPixel[subpixel_counter] << " " << subpixel_counter << " " << cl
0375 << std::endl;
0376 }
0377 }
0378 }
0379
0380 stop = true;
0381 for (unsigned int subcluster_index = 0; subcluster_index < meanExp; subcluster_index++) {
0382 if (std::abs(clx[subcluster_index] - oldclx[subcluster_index]) > 0.01f)
0383 stop = false;
0384 if (std::abs(cly[subcluster_index] - oldcly[subcluster_index]) > 0.01f)
0385 stop = false;
0386 oldclx[subcluster_index] = clx[subcluster_index];
0387 oldcly[subcluster_index] = cly[subcluster_index];
0388 clx[subcluster_index] = 0;
0389 cly[subcluster_index] = 0;
0390 cls[subcluster_index] = 1e-99;
0391 }
0392 for (unsigned int pixel_index = 0; pixel_index < pixels.size(); pixel_index++) {
0393 if (clusterForPixel[pixel_index] < 0)
0394 continue;
0395 if (verbose)
0396 std::cout << "j " << pixel_index << " x " << pixels[pixel_index].second.x << " * y "
0397 << pixels[pixel_index].second.y << " * ADC " << pixels[pixel_index].second.adc << " * W "
0398 << weightOfPixel[pixel_index] << std::endl;
0399 clx[clusterForPixel[pixel_index]] += pixels[pixel_index].second.x * pixels[pixel_index].second.adc;
0400 cly[clusterForPixel[pixel_index]] += pixels[pixel_index].second.y * pixels[pixel_index].second.adc;
0401 cls[clusterForPixel[pixel_index]] += pixels[pixel_index].second.adc;
0402 }
0403 for (unsigned int subcluster_index = 0; subcluster_index < meanExp; subcluster_index++) {
0404 if (cls[subcluster_index] != 0) {
0405 clx[subcluster_index] /= cls[subcluster_index];
0406 cly[subcluster_index] /= cls[subcluster_index];
0407 }
0408 if (verbose)
0409 std::cout << "Center for cluster " << subcluster_index << " x,y " << clx[subcluster_index] << " "
0410 << cly[subcluster_index] << std::endl;
0411 cls[subcluster_index] = 0;
0412 }
0413 }
0414 if (verbose)
0415 std::cout << "maxstep " << remainingSteps << std::endl;
0416
0417 std::vector<std::vector<SiPixelCluster::Pixel>> pixelsForCl(meanExp);
0418 for (int cl = 0; cl < (int)meanExp; cl++) {
0419 for (unsigned int j = 0; j < pixels.size(); j++) {
0420 if (clusterForPixel[j] == cl and pixels[j].second.adc != 0) {
0421
0422
0423
0424 for (unsigned int k = j + 1; k < pixels.size(); k++) {
0425 if (pixels[k].second.adc != 0 and pixels[k].second.x == pixels[j].second.x and
0426 pixels[k].second.y == pixels[j].second.y and clusterForPixel[k] == cl) {
0427 if (verbose)
0428 std::cout << "Resetting all sub-pixel for location " << pixels[k].second.x << ", " << pixels[k].second.y
0429 << " at index " << k << " associated to cl " << clusterForPixel[k] << std::endl;
0430 pixels[j].second.adc += pixels[k].second.adc;
0431 pixels[k].second.adc = 0;
0432 }
0433 }
0434 for (unsigned int p = 0; p < pixels.size(); ++p)
0435 if (verbose)
0436 std::cout << "index, x, y, ADC: " << p << ", " << pixels[p].second.x << ", " << pixels[p].second.y << ", "
0437 << pixels[p].second.adc << " associated to cl " << clusterForPixel[p] << std::endl
0438 << "Adding pixel " << pixels[j].second.x << ", " << pixels[j].second.y << " to cluster " << cl
0439 << std::endl;
0440 pixelsForCl[cl].push_back(pixels[j].second);
0441 }
0442 }
0443 }
0444
0445
0446
0447
0448
0449 for (int cl = 0; cl < (int)meanExp; cl++) {
0450 if (verbose)
0451 std::cout << "Pixels of cl " << cl << " ";
0452 for (unsigned int j = 0; j < pixelsForCl[cl].size(); j++) {
0453 SiPixelCluster::PixelPos newpix(pixelsForCl[cl][j].x, pixelsForCl[cl][j].y);
0454 if (verbose)
0455 std::cout << pixelsForCl[cl][j].x << "," << pixelsForCl[cl][j].y << "|";
0456 if (j == 0) {
0457 output.emplace_back(newpix, pixelsForCl[cl][j].adc);
0458 } else {
0459 output.back().add(newpix, pixelsForCl[cl][j].adc);
0460 }
0461 }
0462 if (verbose)
0463 std::cout << std::endl;
0464 if (!pixelsForCl[cl].empty()) {
0465 if (forceXError_ > 0)
0466 output.back().setSplitClusterErrorX(forceXError_);
0467 if (forceYError_ > 0)
0468 output.back().setSplitClusterErrorY(forceYError_);
0469 }
0470 }
0471
0472
0473
0474
0475
0476 return output;
0477 }
0478
0479 #include "FWCore/PluginManager/interface/ModuleDef.h"
0480 #include "FWCore/Framework/interface/MakerMacros.h"
0481 DEFINE_FWK_MODULE(JetCoreClusterSplitter);