Back to home page

Project CMSSW displayed by LXR
 
 

    

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             // check if the cluster has to be splitted
0128 
0129             bool isEndCap = (std::abs(cPos.z()) > 30.f);  // FIXME: check detID instead!
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             }  // in endcap col/rows are switched
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           // blowup the error if we failed to split a splittable cluster (does
0165           // it ever happen)
0166           c.setSplitClusterErrorX(c.sizeX() * (100.f / 3.f));  // this is not really blowing up .. TODO: tune
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     }  // loop over clusters
0175     std::sort_heap(filler.begin(), filler.end(), [](SiPixelCluster const& cl1, SiPixelCluster const& cl2) {
0176       return cl1.minPixelRow() < cl2.minPixelRow();
0177     });
0178 
0179   }  // loop over det
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   // This function should test several configuration of splitting, and then
0190   // return the one with best chi2
0191 
0192   std::vector<SiPixelCluster> sp = fittingSplit(
0193       aCluster, expectedADC, sizeY, sizeX, jetZOverRho, std::floor(aCluster.charge() / expectedADC + 0.5f));
0194 
0195   // for the config with best chi2
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 // order with fast algo and pick first and second instead?
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   // initial values
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     // Compute all distances
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         //              float sizeX=2;
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     // Compute scores for sequential addition. The first index is the
0331     // distance, in whatever metrics we use, while the second is the
0332     // pixel index w.r.t which the distance is computed.
0333     std::multimap<float, int> scores;
0334 
0335     // Using different rankings to improve convergence (as Giulio proposed)
0336     scores = secondDistScore(distanceMap);
0337 
0338     // Iterate starting from the ones with furthest second best clusters, i.e.
0339     // easy choices
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       // find cluster that is both close and has some charge still to assign
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_);        // 20% uncertainty? realistic from Landau?
0358             float clQest = 1.f / (1.f + std::exp(nsig)) + 1e-6f;  // 1./(1.+exp(x*x-3*3))
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     // Recompute cluster centers
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;  // still moving
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   // accumulate pixel with same cl
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) {  // for each pixel of cluster
0421                                                                      // cl find the other pixels
0422                                                                      // with same x,y and
0423                                                                      // accumulate+reset their adc
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   //    std::vector<std::vector<std::vector<SiPixelCluster::PixelPos *> > >
0446   //pixelMap(meanExp,std::vector<std::vector<SiPixelCluster::PixelPos *>
0447   //>(512,std::vector<SiPixelCluster::Pixel *>(512,0)));
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   //    if(verbose) std::cout << "Weights" << std::endl;
0472   //    if(verbose) print(theWeights,aCluster,1);
0473   //    if(verbose) std::cout << "Unused charge" << std::endl;
0474   //    if(verbose) print(theBufferResidual,aCluster);
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);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!