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