SubCollectionProducers/src/JetCoreClusterSplitter.cc

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 /* px barrel */ && 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             // check if the cluster has to be splitted
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           // blowup the error if we failed to split a splittable cluster (does
0216           // it ever happen)
0217           const float fromCentiToMicro = 1e4;
0218           c.setSplitClusterErrorX(c.sizeX() *
0219                                   (pitchX * fromCentiToMicro / 3.f));  // this is not really blowing up .. TODO: tune
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     }  // loop over clusters
0228     std::sort_heap(filler.begin(), filler.end(), [](SiPixelCluster const& cl1, SiPixelCluster const& cl2) {
0229       return cl1.minPixelRow() < cl2.minPixelRow();
0230     });
0231
0232   }  // loop over det
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   // This function should test several configuration of splitting, and then
0243   // return the one with best chi2
0244
0245   std::vector<SiPixelCluster> sp = fittingSplit(
0246       aCluster, expectedADC, sizeY, sizeX, jetZOverRho, std::floor(aCluster.charge() / expectedADC + 0.5f));
0247
0248   // for the config with best chi2
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 // order with fast algo and pick first and second instead?
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   // initial values
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     // Compute all distances
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         //              float sizeX=2;
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     // Compute scores for sequential addition. The first index is the
0384     // distance, in whatever metrics we use, while the second is the
0385     // pixel index w.r.t which the distance is computed.
0386     std::multimap<float, int> scores;
0387
0388     // Using different rankings to improve convergence (as Giulio proposed)
0389     scores = secondDistScore(distanceMap);
0390
0391     // Iterate starting from the ones with furthest second best clusters, i.e.
0392     // easy choices
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       // find cluster that is both close and has some charge still to assign
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_);        // 20% uncertainty? realistic from Landau?
0411             float clQest = 1.f / (1.f + std::exp(nsig)) + 1e-6f;  // 1./(1.+exp(x*x-3*3))
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     // Recompute cluster centers
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;  // still moving
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   // accumulate pixel with same cl
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) {  // for each pixel of cluster
0474                                                                      // cl find the other pixels
0475                                                                      // with same x,y and
0476                                                                      // accumulate+reset their adc
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   //    std::vector<std::vector<std::vector<SiPixelCluster::PixelPos *> > >
0499   //pixelMap(meanExp,std::vector<std::vector<SiPixelCluster::PixelPos *>
0500   //>(512,std::vector<SiPixelCluster::Pixel *>(512,0)));
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   //    if(verbose) std::cout << "Weights" << std::endl;
0525   //    if(verbose) print(theWeights,aCluster,1);
0526   //    if(verbose) std::cout << "Unused charge" << std::endl;
0527   //    if(verbose) print(theBufferResidual,aCluster);
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);