CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Member Functions | Private Member Functions | Private Attributes
JetCoreClusterSplitter Class Reference
Inheritance diagram for JetCoreClusterSplitter:
edm::stream::EDProducer<> edm::stream::EDProducerBase edm::ProducerBase edm::EDConsumerBase edm::ProductRegistryHelper

Public Member Functions

 JetCoreClusterSplitter (const edm::ParameterSet &iConfig)
 
void produce (edm::Event &iEvent, const edm::EventSetup &iSetup)
 
 ~JetCoreClusterSplitter ()
 
- Public Member Functions inherited from edm::stream::EDProducer<>
 EDProducer ()=default
 
- Public Member Functions inherited from edm::stream::EDProducerBase
 EDProducerBase ()
 
ModuleDescription const & moduleDescription () const
 
virtual ~EDProducerBase ()
 
- Public Member Functions inherited from edm::ProducerBase
void callWhenNewProductsRegistered (std::function< void(BranchDescription const &)> const &func)
 
std::vector
< edm::ProductResolverIndex >
const & 
indiciesForPutProducts (BranchType iBranchType) const
 
 ProducerBase ()
 
void registerProducts (ProducerBase *, ProductRegistry *, ModuleDescription const &)
 
std::function< void(BranchDescription
const &)> 
registrationCallback () const
 used by the fwk to register list of products More...
 
void resolvePutIndicies (BranchType iBranchType, std::unordered_multimap< std::string, edm::ProductResolverIndex > const &iIndicies, std::string const &moduleLabel)
 
virtual ~ProducerBase () noexcept(false)
 
- Public Member Functions inherited from edm::EDConsumerBase
std::vector< ConsumesInfoconsumesInfo () const
 
 EDConsumerBase ()
 
 EDConsumerBase (EDConsumerBase const &)=delete
 
 EDConsumerBase (EDConsumerBase &&)=default
 
ProductResolverIndexAndSkipBit indexFrom (EDGetToken, BranchType, TypeID const &) const
 
void itemsMayGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const
 
void itemsToGet (BranchType, std::vector< ProductResolverIndexAndSkipBit > &) const
 
std::vector
< ProductResolverIndexAndSkipBit >
const & 
itemsToGetFromEvent () const
 
void labelsForToken (EDGetToken iToken, Labels &oLabels) const
 
void modulesWhoseProductsAreConsumed (std::vector< ModuleDescription const * > &modules, ProductRegistry const &preg, std::map< std::string, ModuleDescription const * > const &labelsToDesc, std::string const &processName) const
 
EDConsumerBase const & operator= (EDConsumerBase const &)=delete
 
EDConsumerBaseoperator= (EDConsumerBase &&)=default
 
bool registeredToConsume (ProductResolverIndex, bool, BranchType) const
 
bool registeredToConsumeMany (TypeID const &, BranchType) const
 
void updateLookup (BranchType iBranchType, ProductResolverIndexHelper const &, bool iPrefetchMayGet)
 
virtual ~EDConsumerBase () noexcept(false)
 

Private Member Functions

std::pair< float, float > closestClusters (const std::vector< float > &distanceMap)
 
std::multimap< float, int > distScore (const std::vector< std::vector< float > > &distanceMap)
 
std::vector< SiPixelClusterfittingSplit (const SiPixelCluster &aCluster, float expectedADC, int sizeY, int sizeX, float jetZOverRho, unsigned int nSplitted)
 
std::multimap< float, int > secondDistDiffScore (const std::vector< std::vector< float > > &distanceMap)
 
std::multimap< float, int > secondDistScore (const std::vector< std::vector< float > > &distanceMap)
 
bool split (const SiPixelCluster &aCluster, edmNew::DetSetVector< SiPixelCluster >::FastFiller &filler, float expectedADC, int sizeY, int sizeX, float jetZOverRho)
 

Private Attributes

double centralMIPCharge_
 
double chargeFracMin_
 
double chargePerUnit_
 
edm::EDGetTokenT< edm::View
< reco::Candidate > > 
cores_
 
double deltaR_
 
double forceXError_
 
double forceYError_
 
double fractionalWidth_
 
edm::EDGetTokenT
< edmNew::DetSetVector
< SiPixelCluster > > 
pixelClusters_
 
std::string pixelCPE_
 
double ptMin_
 
bool verbose
 
edm::EDGetTokenT
< reco::VertexCollection
vertices_
 

Additional Inherited Members

- Public Types inherited from edm::stream::EDProducer<>
typedef CacheContexts< T...> CacheTypes
 
typedef CacheTypes::GlobalCache GlobalCache
 
typedef AbilityChecker< T...> HasAbility
 
typedef
CacheTypes::LuminosityBlockCache 
LuminosityBlockCache
 
typedef
LuminosityBlockContextT
< LuminosityBlockCache,
RunCache, GlobalCache
LuminosityBlockContext
 
typedef
CacheTypes::LuminosityBlockSummaryCache 
LuminosityBlockSummaryCache
 
typedef CacheTypes::RunCache RunCache
 
typedef RunContextT< RunCache,
GlobalCache
RunContext
 
typedef CacheTypes::RunSummaryCache RunSummaryCache
 
- Public Types inherited from edm::stream::EDProducerBase
typedef EDProducerAdaptorBase ModuleType
 
- Public Types inherited from edm::ProducerBase
typedef
ProductRegistryHelper::TypeLabelList 
TypeLabelList
 
- Public Types inherited from edm::EDConsumerBase
typedef ProductLabels Labels
 
- Static Public Member Functions inherited from edm::stream::EDProducerBase
static const std::string & baseType ()
 
static void fillDescriptions (ConfigurationDescriptions &descriptions)
 
static void prevalidate (ConfigurationDescriptions &descriptions)
 
- Protected Member Functions inherited from edm::EDConsumerBase
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType > consumes (edm::InputTag const &tag)
 
EDGetToken consumes (const TypeToGet &id, edm::InputTag const &tag)
 
template<BranchType B>
EDGetToken consumes (TypeToGet const &id, edm::InputTag const &tag)
 
ConsumesCollector consumesCollector ()
 Use a ConsumesCollector to gather consumes information from helper functions. More...
 
template<typename ProductType , BranchType B = InEvent>
void consumesMany ()
 
void consumesMany (const TypeToGet &id)
 
template<BranchType B>
void consumesMany (const TypeToGet &id)
 
template<typename ProductType , BranchType B = InEvent>
EDGetTokenT< ProductType > mayConsume (edm::InputTag const &tag)
 
EDGetToken mayConsume (const TypeToGet &id, edm::InputTag const &tag)
 
template<BranchType B>
EDGetToken mayConsume (const TypeToGet &id, edm::InputTag const &tag)
 

Detailed Description

Definition at line 26 of file JetCoreClusterSplitter.cc.

Constructor & Destructor Documentation

JetCoreClusterSplitter::JetCoreClusterSplitter ( const edm::ParameterSet iConfig)

Definition at line 63 of file JetCoreClusterSplitter.cc.

64  : verbose(iConfig.getParameter<bool>("verbose")),
65  pixelCPE_(iConfig.getParameter<std::string>("pixelCPE")),
66  ptMin_(iConfig.getParameter<double>("ptMin")),
67  deltaR_(iConfig.getParameter<double>("deltaRmax")),
68  chargeFracMin_(iConfig.getParameter<double>("chargeFractionMin")),
70  iConfig.getParameter<edm::InputTag>("pixelClusters"))),
71  vertices_(consumes<reco::VertexCollection>(
72  iConfig.getParameter<edm::InputTag>("vertices"))),
74  iConfig.getParameter<edm::InputTag>("cores"))),
75  forceXError_(iConfig.getParameter<double>("forceXError")),
76  forceYError_(iConfig.getParameter<double>("forceYError")),
77  fractionalWidth_(iConfig.getParameter<double>("fractionalWidth")),
78  chargePerUnit_(iConfig.getParameter<double>("chargePerUnit")),
79  centralMIPCharge_(iConfig.getParameter<double>("centralMIPCharge"))
80 
81 {
82  produces<edmNew::DetSetVector<SiPixelCluster> >();
83 }
T getParameter(std::string const &) const
edm::EDGetTokenT< edmNew::DetSetVector< SiPixelCluster > > pixelClusters_
edm::EDGetTokenT< reco::VertexCollection > vertices_
edm::EDGetTokenT< edm::View< reco::Candidate > > cores_
EDGetTokenT< ProductType > consumes(edm::InputTag const &tag)
JetCoreClusterSplitter::~JetCoreClusterSplitter ( )

Definition at line 85 of file JetCoreClusterSplitter.cc.

85 {}

Member Function Documentation

std::pair< float, float > JetCoreClusterSplitter::closestClusters ( const std::vector< float > &  distanceMap)
private

Definition at line 216 of file JetCoreClusterSplitter.cc.

References i.

Referenced by distScore(), secondDistDiffScore(), and secondDistScore().

217  {
218  float minDist = 9e99;
219  float secondMinDist = 9e99;
220  for (unsigned int i = 0; i < distanceMap.size(); i++) {
221  float dist = distanceMap[i];
222  if (dist < minDist) {
223  secondMinDist = minDist;
224  minDist = dist;
225  } else if (dist < secondMinDist) {
226  secondMinDist = dist;
227  }
228  }
229  return std::pair<float, float>(minDist, secondMinDist);
230 }
int i
Definition: DBlmapReader.cc:9
std::multimap< float, int > JetCoreClusterSplitter::distScore ( const std::vector< std::vector< float > > &  distanceMap)
private

Definition at line 252 of file JetCoreClusterSplitter.cc.

References closestClusters(), ztail::d, and j.

253  {
254  std::multimap<float, int> scores;
255  for (unsigned int j = 0; j < distanceMap.size(); j++) {
256  std::pair<float, float> d = closestClusters(distanceMap[j]);
257  scores.insert(std::pair<float, int>(-d.first, j));
258  }
259  return scores;
260 }
tuple d
Definition: ztail.py:151
int j
Definition: DBlmapReader.cc:9
std::pair< float, float > closestClusters(const std::vector< float > &distanceMap)
std::vector< SiPixelCluster > JetCoreClusterSplitter::fittingSplit ( const SiPixelCluster aCluster,
float  expectedADC,
int  sizeY,
int  sizeX,
float  jetZOverRho,
unsigned int  nSplitted 
)
private

Definition at line 262 of file JetCoreClusterSplitter.cc.

References funct::abs(), ecalMGPA::adc(), centralMIPCharge_, chargePerUnit_, GetRecoTauVFromDQM_MC_cff::cl, gather_cfg::cout, alignCSCRings::e, create_public_lumi_plots::exp, forceXError_, forceYError_, fractionalWidth_, i, j, relval_2017::k, convertSQLitetoXML_cfg::output, AlCaHLTBitMon_ParallelJobs::p, SiPixelCluster::pixels(), edm::second(), secondDistScore(), findQualityFiles::size, mathSSE::sqrt(), x, and y.

Referenced by split().

264  {
265  std::vector<SiPixelCluster> output;
266 
267  unsigned int meanExp = nSplitted;
268  if (meanExp <= 1) {
269  output.push_back(aCluster);
270  return output;
271  }
272 
273  std::vector<float> clx(meanExp);
274  std::vector<float> cly(meanExp);
275  std::vector<float> cls(meanExp);
276  std::vector<float> oldclx(meanExp);
277  std::vector<float> oldcly(meanExp);
278  std::vector<SiPixelCluster::Pixel> originalpixels = aCluster.pixels();
279  std::vector<std::pair<int, SiPixelCluster::Pixel> > pixels;
280  for (unsigned int j = 0; j < originalpixels.size(); j++) {
281  int sub = originalpixels[j].adc / chargePerUnit_ * expectedADC /
283  if (sub < 1) sub = 1;
284  int perDiv = originalpixels[j].adc / sub;
285  if (verbose)
286  std::cout << "Splitting " << j << " in [ " << pixels.size() << " , "
287  << pixels.size() + sub << " ], expected numb of clusters: "
288  << meanExp << " original pixel (x,y) "
289  << originalpixels[j].x << " " << originalpixels[j].y
290  << " sub " << sub << std::endl;
291  for (int k = 0; k < sub; k++) {
292  if (k == sub - 1) perDiv = originalpixels[j].adc - perDiv * k;
293  pixels.push_back(std::make_pair(j, SiPixelCluster::Pixel(originalpixels[j].x,
294  originalpixels[j].y, perDiv)));
295  }
296  }
297  std::vector<int> clusterForPixel(pixels.size());
298  // initial values
299  for (unsigned int j = 0; j < meanExp; j++) {
300  oldclx[j] = -999;
301  oldcly[j] = -999;
302  clx[j] = originalpixels[0].x + j;
303  cly[j] = originalpixels[0].y + j;
304  cls[j] = 0;
305  }
306  bool stop = false;
307  int remainingSteps = 100;
308  while (!stop && remainingSteps > 0) {
309  remainingSteps--;
310  // Compute all distances
311  std::vector<std::vector<float> > distanceMapX(originalpixels.size(),
312  std::vector<float>(meanExp));
313  std::vector<std::vector<float> > distanceMapY(originalpixels.size(),
314  std::vector<float>(meanExp));
315  std::vector<std::vector<float> > distanceMap(originalpixels.size(),
316  std::vector<float>(meanExp));
317  for (unsigned int j = 0; j < originalpixels.size(); j++) {
318  if (verbose)
319  std::cout << "Original Pixel pos " << j << " " << pixels[j].second.x << " "
320  << pixels[j].second.y << std::endl;
321  for (unsigned int i = 0; i < meanExp; i++) {
322  distanceMapX[j][i] = 1. * originalpixels[j].x - clx[i];
323  distanceMapY[j][i] = 1. * originalpixels[j].y - cly[i];
324  float dist = 0;
325  // float sizeX=2;
326  if (std::abs(distanceMapX[j][i]) > sizeX / 2.) {
327  dist += (std::abs(distanceMapX[j][i]) - sizeX / 2. + 1) *
328  (std::abs(distanceMapX[j][i]) - sizeX / 2. + 1);
329  } else {
330  dist +=
331  distanceMapX[j][i] / sizeX * 2 * distanceMapX[j][i] / sizeX * 2;
332  }
333 
334  if (std::abs(distanceMapY[j][i]) > sizeY / 2.) {
335  dist += 1. * (std::abs(distanceMapY[j][i]) - sizeY / 2. + 1.) *
336  (std::abs(distanceMapY[j][i]) - sizeY / 2. + 1.);
337  } else {
338  dist += 1. * distanceMapY[j][i] / sizeY * 2. * distanceMapY[j][i] /
339  sizeY * 2.;
340  }
341  distanceMap[j][i] = sqrt(dist);
342  if (verbose)
343  std::cout << "Cluster " << i << " Original Pixel " << j
344  << " distances: " << distanceMapX[j][i] << " "
345  << distanceMapY[j][i] << " " << distanceMap[j][i]
346  << std::endl;
347  }
348  }
349  // Compute scores for sequential addition. The first index is the
350  // distance, in whatever metrics we use, while the second is the
351  // pixel index w.r.t which the distance is computed.
352  std::multimap<float, int> scores;
353 
354  // Using different rankings to improve convergence (as Giulio proposed)
355  scores = secondDistScore(distanceMap);
356 
357  // Iterate starting from the ones with furthest second best clusters, i.e.
358  // easy choices
359  std::vector<float> weightOfPixel(pixels.size());
360  for (std::multimap<float, int>::iterator it = scores.begin();
361  it != scores.end(); it++) {
362  int pixel_index = it->second;
363  if (verbose)
364  std::cout << "Original Pixel " << pixel_index << " with score " << it->first << std::endl;
365  // find cluster that is both close and has some charge still to assign
366  int subpixel_counter = 0;
367  for (auto subpixel = pixels.begin(); subpixel != pixels.end();
368  ++subpixel, ++subpixel_counter) {
369  if (subpixel->first > pixel_index) {
370  break;
371  } else if (subpixel->first != pixel_index) {
372  continue;
373  } else {
374  float maxEst = 0;
375  int cl = -1;
376  for (unsigned int subcluster_index = 0;
377  subcluster_index < meanExp; subcluster_index++) {
378  float nsig =
379  (cls[subcluster_index] - expectedADC) /
380  (expectedADC *
381  fractionalWidth_); // 20% uncertainty? realistic from Landau?
382  float clQest = 1. / (1. + exp(nsig)) + 1e-6; // 1./(1.+exp(x*x-3*3))
383  float clDest = 1. / (distanceMap[pixel_index][subcluster_index] + 0.05);
384 
385  if (verbose)
386  std::cout << " Q: " << clQest << " D: " << clDest << " "
387  << distanceMap[pixel_index][subcluster_index] << std::endl;
388  float est = clQest * clDest;
389  if (est > maxEst) {
390  cl = subcluster_index;
391  maxEst = est;
392  }
393  }
394  cls[cl] += subpixel->second.adc;
395  clusterForPixel[subpixel_counter] = cl;
396  weightOfPixel[subpixel_counter] = maxEst;
397  if (verbose)
398  std::cout << "Pixel weight j cl " << weightOfPixel[subpixel_counter]
399  << " " << subpixel_counter
400  << " " << cl << std::endl;
401  }
402  }
403  }
404  // Recompute cluster centers
405  stop = true;
406  for (unsigned int subcluster_index = 0;
407  subcluster_index < meanExp; subcluster_index++) {
408  if (std::abs(clx[subcluster_index] - oldclx[subcluster_index]) > 0.01)
409  stop = false; // still moving
410  if (std::abs(cly[subcluster_index] - oldcly[subcluster_index]) > 0.01)
411  stop = false;
412  oldclx[subcluster_index] = clx[subcluster_index];
413  oldcly[subcluster_index] = cly[subcluster_index];
414  clx[subcluster_index] = 0;
415  cly[subcluster_index] = 0;
416  cls[subcluster_index] = 1e-99;
417  }
418  for (unsigned int pixel_index = 0;
419  pixel_index < pixels.size(); pixel_index++) {
420  if (clusterForPixel[pixel_index] < 0) continue;
421  if (verbose)
422  std::cout << "j " << pixel_index << " x " << pixels[pixel_index].second.x << " * y "
423  << pixels[pixel_index].second.y << " * ADC "
424  << pixels[pixel_index].second.adc << " * W "
425  << weightOfPixel[pixel_index] << std::endl;
426  clx[clusterForPixel[pixel_index]] += pixels[pixel_index].second.x * pixels[pixel_index].second.adc;
427  cly[clusterForPixel[pixel_index]] += pixels[pixel_index].second.y * pixels[pixel_index].second.adc;
428  cls[clusterForPixel[pixel_index]] += pixels[pixel_index].second.adc;
429  }
430  for (unsigned int subcluster_index = 0;
431  subcluster_index < meanExp; subcluster_index++) {
432  if (cls[subcluster_index] != 0) {
433  clx[subcluster_index] /= cls[subcluster_index];
434  cly[subcluster_index] /= cls[subcluster_index];
435  }
436  if (verbose)
437  std::cout << "Center for cluster " << subcluster_index << " x,y "
438  << clx[subcluster_index] << " "
439  << cly[subcluster_index] << std::endl;
440  cls[subcluster_index] = 0;
441  }
442  }
443  if (verbose) std::cout << "maxstep " << remainingSteps << std::endl;
444  // accumulate pixel with same cl
445  std::vector<std::vector<SiPixelCluster::Pixel> > pixelsForCl(meanExp);
446  for (int cl = 0; cl < (int)meanExp; cl++) {
447  for (unsigned int j = 0; j < pixels.size(); j++) {
448  if (clusterForPixel[j] == cl and
449  pixels[j].second.adc != 0) { // for each pixel of cluster
450  // cl find the other pixels
451  // with same x,y and
452  // accumulate+reset their adc
453  for (unsigned int k = j + 1; k < pixels.size(); k++) {
454  if (pixels[k].second.adc != 0
455  and pixels[k].second.x == pixels[j].second.x
456  and pixels[k].second.y == pixels[j].second.y
457  and clusterForPixel[k] == cl) {
458  if (verbose)
459  std::cout << "Resetting all sub-pixel for location "
460  << pixels[k].second.x << ", " << pixels[k].second.y
461  << " at index " << k << " associated to cl "
462  << clusterForPixel[k] << std::endl;
463  pixels[j].second.adc += pixels[k].second.adc;
464  pixels[k].second.adc = 0;
465  }
466  }
467  for (unsigned int p = 0; p < pixels.size(); ++p)
468  if (verbose)
469  std::cout << "index, x, y, ADC: " << p << ", "
470  << pixels[p].second.x << ", " << pixels[p].second.y
471  << ", " << pixels[p].second.adc
472  << " associated to cl " << clusterForPixel[p] << std::endl
473  << "Adding pixel " << pixels[j].second.x << ", " << pixels[j].second.y
474  << " to cluster " << cl << std::endl;
475  pixelsForCl[cl].push_back(pixels[j].second);
476  }
477  }
478  }
479 
480  // std::vector<std::vector<std::vector<SiPixelCluster::PixelPos *> > >
481  //pixelMap(meanExp,std::vector<std::vector<SiPixelCluster::PixelPos *>
482  //>(512,std::vector<SiPixelCluster::Pixel *>(512,0)));
483 
484  for (int cl = 0; cl < (int)meanExp; cl++) {
485  if (verbose) std::cout << "Pixels of cl " << cl << " ";
486  for (unsigned int j = 0; j < pixelsForCl[cl].size(); j++) {
487  SiPixelCluster::PixelPos newpix(pixelsForCl[cl][j].x,
488  pixelsForCl[cl][j].y);
489  if (verbose)
490  std::cout << pixelsForCl[cl][j].x << "," << pixelsForCl[cl][j].y << "|";
491  if (j==0) {
492  output.emplace_back(newpix, pixelsForCl[cl][j].adc);
493  } else {
494  output.back().add(newpix, pixelsForCl[cl][j].adc);
495  }
496  }
497  if (verbose) std::cout << std::endl;
498  if (pixelsForCl[cl].size() > 0) {
499  if (forceXError_ > 0) output.back().setSplitClusterErrorX(forceXError_);
500  if (forceYError_ > 0) output.back().setSplitClusterErrorY(forceYError_);
501  }
502  }
503  // if(verbose) std::cout << "Weights" << std::endl;
504  // if(verbose) print(theWeights,aCluster,1);
505  // if(verbose) std::cout << "Unused charge" << std::endl;
506  // if(verbose) print(theBufferResidual,aCluster);
507 
508  return output;
509 }
int adc(sample_type sample)
get the ADC sample (12 bits)
int i
Definition: DBlmapReader.cc:9
U second(std::pair< T, U > const &p)
T sqrt(T t)
Definition: SSEVec.h:18
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int j
Definition: DBlmapReader.cc:9
tuple cout
Definition: gather_cfg.py:145
std::multimap< float, int > secondDistScore(const std::vector< std::vector< float > > &distanceMap)
tuple size
Write out results.
const std::vector< Pixel > pixels() const
void JetCoreClusterSplitter::produce ( edm::Event iEvent,
const edm::EventSetup iSetup 
)
virtual

Implements edm::stream::EDProducerBase.

Definition at line 92 of file JetCoreClusterSplitter.cc.

References edmNew::DetSetVector< T >::begin(), EnergyCorrector::c, centralMIPCharge_, SiPixelCluster::charge(), chargeFracMin_, HLT_FULL_cff::cores, cores_, gather_cfg::cout, deltaR(), deltaR_, edmNew::DetSetVector< T >::end(), plotBeamSpotDB::first, geometry, edm::EventSetup::get(), edm::Event::getByToken(), edmNew::DetSetVector< T >::id(), metsig::jet, reco::Candidate::momentum(), eostools::move(), convertSQLitetoXML_cfg::output, pixelClusters_, pixelCPE_, createTree::pp, ptMin_, edm::Event::put(), MetAnalyzer::pv(), reco::Candidate::px(), reco::Candidate::py(), reco::Candidate::pz(), SiPixelCluster::setSplitClusterErrorX(), SiPixelCluster::setSplitClusterErrorY(), SiPixelCluster::sizeX(), SiPixelCluster::sizeY(), split(), mathSSE::sqrt(), GeomDet::surface(), Surface::toGlobal(), HLT_FULL_cff::vertices, vertices_, and PV3DBase< T, PVType, FrameType >::z().

93  {
94  using namespace edm;
97 
99  iEvent.getByToken(pixelClusters_, inputPixelClusters);
100 
102  iEvent.getByToken(vertices_, vertices);
103  const reco::Vertex& pv = (*vertices)[0];
104 
106  iEvent.getByToken(cores_, cores);
107 
110  iSetup.get<TkPixelCPERecord>().get(pixelCPE_, pe);
111  pp = pe.product();
112 
113  auto output = std::make_unique<edmNew::DetSetVector<SiPixelCluster>>();
114 
116  inputPixelClusters->begin();
117  for (; detIt != inputPixelClusters->end(); detIt++) {
119  detIt->id());
120  const edmNew::DetSet<SiPixelCluster>& detset = *detIt;
121  const GeomDet* det = geometry->idToDet(detset.id());
123  detset.begin();
124  cluster != detset.end(); cluster++) {
125  const SiPixelCluster& aCluster = *cluster;
126  bool hasBeenSplit = false;
127  bool shouldBeSplit = false;
128  GlobalPoint cPos = det->surface().toGlobal(
129  pp->localParametersV(aCluster,
130  (*geometry->idToDetUnit(detIt->id())))[0].first);
131  GlobalPoint ppv(pv.position().x(), pv.position().y(), pv.position().z());
132  GlobalVector clusterDir = cPos - ppv;
133  for (unsigned int ji = 0; ji < cores->size(); ji++) {
134  if ((*cores)[ji].pt() > ptMin_) {
135  const reco::Candidate& jet = (*cores)[ji];
136  GlobalVector jetDir(jet.px(), jet.py(), jet.pz());
137  if (Geom::deltaR(jetDir, clusterDir) < deltaR_) {
138  // check if the cluster has to be splitted
139 
140  bool isEndCap =
141  (fabs(cPos.z()) > 30); // FIXME: check detID instead!
142  float jetZOverRho = jet.momentum().Z() / jet.momentum().Rho();
143  if (isEndCap)
144  jetZOverRho = jet.momentum().Rho() / jet.momentum().Z();
145  float expSizeY =
146  fabs(sqrt(1.3 * 1.3 + 1.9 * 1.9 * jetZOverRho * jetZOverRho));
147  if (expSizeY < 1) expSizeY = 1.;
148  float expSizeX = 1.5;
149  if (isEndCap) {
150  expSizeX = expSizeY;
151  expSizeY = 1.5;
152  } // in endcap col/rows are switched
153  float expCharge =
154  sqrt(1.08 + jetZOverRho * jetZOverRho) * centralMIPCharge_;
155 
156  if (aCluster.charge() > expCharge * chargeFracMin_ &&
157  (aCluster.sizeX() > expSizeX + 1 ||
158  aCluster.sizeY() > expSizeY + 1)) {
159  shouldBeSplit = true;
160  if (verbose)
161  std::cout << "Trying to split: charge and deltaR "
162  << aCluster.charge() << " "
163  << Geom::deltaR(jetDir, clusterDir) << " size x y "
164  << aCluster.sizeX() << " " << aCluster.sizeY()
165  << " exp. size (x,y) "
166  << expSizeX << " " << expSizeY
167  << " detid " << detIt->id() << std::endl;
168  if (verbose)
169  std::cout << "jetZOverRho=" << jetZOverRho << std::endl;
170 
171  if (split(aCluster, filler, expCharge, expSizeY, expSizeX,
172  jetZOverRho)) {
173  hasBeenSplit = true;
174  }
175  }
176  }
177  }
178  }
179  if (!hasBeenSplit) {
180  SiPixelCluster c = aCluster;
181  if (shouldBeSplit) {
182  // blowup the error if we failed to split a splittable cluster (does
183  // it ever happen)
185  c.sizeX() * 100. /
186  3.); // this is not really blowing up .. TODO: tune
187  c.setSplitClusterErrorY(c.sizeY() * 150. / 3.);
188  }
189  filler.push_back(c);
190  }
191  }
192  }
193  iEvent.put(std::move(output));
194 }
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:106
float charge() const
boost::transform_iterator< IterHelp, const_IdIter > const_iterator
const_iterator end(bool update=false) const
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:122
void setSplitClusterErrorY(float erry)
tuple pp
Definition: createTree.py:15
edm::EDGetTokenT< edmNew::DetSetVector< SiPixelCluster > > pixelClusters_
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:457
virtual double pz() const =0
z coordinate of momentum vector
bool split(const SiPixelCluster &aCluster, edmNew::DetSetVector< SiPixelCluster >::FastFiller &filler, float expectedADC, int sizeY, int sizeX, float jetZOverRho)
data_type const * const_iterator
Definition: DetSetNew.h:30
edm::EDGetTokenT< reco::VertexCollection > vertices_
id_type id(size_t cell) const
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:42
edm::EDGetTokenT< edm::View< reco::Candidate > > cores_
virtual Vector momentum() const =0
spatial momentum vector
T sqrt(T t)
Definition: SSEVec.h:18
T z() const
Definition: PV3DBase.h:64
def move
Definition: eostools.py:510
virtual double py() const =0
y coordinate of momentum vector
virtual double px() const =0
x coordinate of momentum vector
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
void setSplitClusterErrorX(float errx)
const T & get() const
Definition: EventSetup.h:56
int sizeY() const
ESHandle< TrackerGeometry > geometry
Pixel cluster – collection of neighboring pixels above threshold.
tuple cout
Definition: gather_cfg.py:145
int sizeX() const
const_iterator begin(bool update=false) const
std::multimap< float, int > JetCoreClusterSplitter::secondDistDiffScore ( const std::vector< std::vector< float > > &  distanceMap)
private

Definition at line 232 of file JetCoreClusterSplitter.cc.

References closestClusters(), ztail::d, and j.

233  {
234  std::multimap<float, int> scores;
235  for (unsigned int j = 0; j < distanceMap.size(); j++) {
236  std::pair<float, float> d = closestClusters(distanceMap[j]);
237  scores.insert(std::pair<float, int>(d.second - d.first, j));
238  }
239  return scores;
240 }
tuple d
Definition: ztail.py:151
int j
Definition: DBlmapReader.cc:9
std::pair< float, float > closestClusters(const std::vector< float > &distanceMap)
std::multimap< float, int > JetCoreClusterSplitter::secondDistScore ( const std::vector< std::vector< float > > &  distanceMap)
private

Definition at line 242 of file JetCoreClusterSplitter.cc.

References closestClusters(), ztail::d, and j.

Referenced by fittingSplit().

243  {
244  std::multimap<float, int> scores;
245  for (unsigned int j = 0; j < distanceMap.size(); j++) {
246  std::pair<float, float> d = closestClusters(distanceMap[j]);
247  scores.insert(std::pair<float, int>(-d.second, j));
248  }
249  return scores;
250 }
tuple d
Definition: ztail.py:151
int j
Definition: DBlmapReader.cc:9
std::pair< float, float > closestClusters(const std::vector< float > &distanceMap)
bool JetCoreClusterSplitter::split ( const SiPixelCluster aCluster,
edmNew::DetSetVector< SiPixelCluster >::FastFiller &  filler,
float  expectedADC,
int  sizeY,
int  sizeX,
float  jetZOverRho 
)
private

Definition at line 196 of file JetCoreClusterSplitter.cc.

References SiPixelCluster::charge(), fittingSplit(), i, and edmNew::DetSetVector< T >::FastFiller::push_back().

Referenced by produce().

199  {
200  // This function should test several configuration of splitting, and then
201  // return the one with best chi2
202 
203  std::vector<SiPixelCluster> sp =
204  fittingSplit(aCluster, expectedADC, sizeY, sizeX, jetZOverRho,
205  floor(aCluster.charge() / expectedADC + 0.5));
206 
207  // for the config with best chi2
208  for (unsigned int i = 0; i < sp.size(); i++) {
209  filler.push_back(sp[i]);
210  }
211 
212  return (sp.size() > 0);
213 }
int i
Definition: DBlmapReader.cc:9
float charge() const
void push_back(data_type const &d)
std::vector< SiPixelCluster > fittingSplit(const SiPixelCluster &aCluster, float expectedADC, int sizeY, int sizeX, float jetZOverRho, unsigned int nSplitted)
Definition: sp.h:21

Member Data Documentation

double JetCoreClusterSplitter::centralMIPCharge_
private

Definition at line 60 of file JetCoreClusterSplitter.cc.

Referenced by fittingSplit(), and produce().

double JetCoreClusterSplitter::chargeFracMin_
private

Definition at line 52 of file JetCoreClusterSplitter.cc.

Referenced by produce().

double JetCoreClusterSplitter::chargePerUnit_
private

Definition at line 59 of file JetCoreClusterSplitter.cc.

Referenced by fittingSplit().

edm::EDGetTokenT<edm::View<reco::Candidate> > JetCoreClusterSplitter::cores_
private

Definition at line 55 of file JetCoreClusterSplitter.cc.

Referenced by produce().

double JetCoreClusterSplitter::deltaR_
private

Definition at line 51 of file JetCoreClusterSplitter.cc.

Referenced by produce().

double JetCoreClusterSplitter::forceXError_
private

Definition at line 56 of file JetCoreClusterSplitter.cc.

Referenced by fittingSplit().

double JetCoreClusterSplitter::forceYError_
private

Definition at line 57 of file JetCoreClusterSplitter.cc.

Referenced by fittingSplit().

double JetCoreClusterSplitter::fractionalWidth_
private

Definition at line 58 of file JetCoreClusterSplitter.cc.

Referenced by fittingSplit().

edm::EDGetTokenT<edmNew::DetSetVector<SiPixelCluster> > JetCoreClusterSplitter::pixelClusters_
private

Definition at line 53 of file JetCoreClusterSplitter.cc.

Referenced by produce().

std::string JetCoreClusterSplitter::pixelCPE_
private

Definition at line 49 of file JetCoreClusterSplitter.cc.

Referenced by produce().

double JetCoreClusterSplitter::ptMin_
private

Definition at line 50 of file JetCoreClusterSplitter.cc.

Referenced by produce().

bool JetCoreClusterSplitter::verbose
private

Definition at line 48 of file JetCoreClusterSplitter.cc.

edm::EDGetTokenT<reco::VertexCollection> JetCoreClusterSplitter::vertices_
private

Definition at line 54 of file JetCoreClusterSplitter.cc.

Referenced by produce().