Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-10-23 22:48:11 

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