CMS 3D CMS Logo

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
 
void convertCurrentProcessAlias (std::string const &processName)
 Convert "@currentProcess" in InputTag process names to the actual current process name. More...
 
 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 & itemsToGetFrom (BranchType iType) 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::VertexCollectionvertices_
 

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, GlobalCacheLuminosityBlockContext
 
typedef CacheTypes::LuminosityBlockSummaryCache LuminosityBlockSummaryCache
 
typedef CacheTypes::RunCache RunCache
 
typedef RunContextT< RunCache, GlobalCacheRunContext
 
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 221 of file JetCoreClusterSplitter.cc.

References mps_fire::i, and hpstanc_transforms::max.

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

222  {
223  float minDist = std::numeric_limits<float>::max();
224  float secondMinDist = std::numeric_limits<float>::max();
225  for (unsigned int i = 0; i < distanceMap.size(); i++) {
226  float dist = distanceMap[i];
227  if (dist < minDist) {
228  secondMinDist = minDist;
229  minDist = dist;
230  } else if (dist < secondMinDist) {
231  secondMinDist = dist;
232  }
233  }
234  return std::pair<float, float>(minDist, secondMinDist);
235 }
std::multimap< float, int > JetCoreClusterSplitter::distScore ( const std::vector< std::vector< float > > &  distanceMap)
private

Definition at line 257 of file JetCoreClusterSplitter.cc.

References closestClusters(), and edmIntegrityCheck::d.

258  {
259  std::multimap<float, int> scores;
260  for (unsigned int j = 0; j < distanceMap.size(); j++) {
261  std::pair<float, float> d = closestClusters(distanceMap[j]);
262  scores.insert(std::pair<float, int>(-d.first, j));
263  }
264  return scores;
265 }
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 267 of file JetCoreClusterSplitter.cc.

References funct::abs(), ecalMGPA::adc(), centralMIPCharge_, chargePerUnit_, GetRecoTauVFromDQM_MC_cff::cl, gather_cfg::cout, DEFINE_FWK_MODULE, MillePedeFileConverter_cfg::e, JetChargeProducer_cfi::exp, f, forceXError_, forceYError_, fractionalWidth_, mps_fire::i, createfilelist::int, gen::k, convertSQLitetoXML_cfg::output, AlCaHLTBitMon_ParallelJobs::p, SiPixelCluster::pixels(), edm::second(), secondDistScore(), findQualityFiles::size, RecoTauValidation_cfi::sizeX, RecoTauValidation_cfi::sizeY, mathSSE::sqrt(), x, and y.

Referenced by split().

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

Definition at line 92 of file JetCoreClusterSplitter.cc.

References funct::abs(), edmNew::DetSetVector< T >::begin(), EnergyCorrector::c, centralMIPCharge_, SiPixelCluster::charge(), chargeFracMin_, GetRecoTauVFromDQM_MC_cff::cl2, HIInitialJetCoreClusterSplitting_cff::cores, cores_, gather_cfg::cout, deltaR(), deltaR_, edmNew::DetSetVector< T >::end(), f, objects.autophobj::filler, plotBeamSpotDB::first, geometry, edm::EventSetup::get(), edm::Event::getByToken(), edmNew::DetSetVector< T >::id(), GlobalTrackingGeometry::idToDet(), GlobalTrackingGeometry::idToDetUnit(), metsig::jet, PixelClusterParameterEstimator::localParametersV(), SiPixelCluster::minPixelRow(), reco::Candidate::momentum(), eostools::move(), convertSQLitetoXML_cfg::output, pixelClusters_, pixelCPE_, reco::Vertex::position(), createTree::pp, edm::ESHandle< T >::product(), 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(), ecalDrivenElectronSeedsParameters_cff::vertices, vertices_, and PV3DBase< T, PVType, FrameType >::z().

93  {
94  using namespace edm;
96  iSetup.get<GlobalTrackingGeometryRecord>().get(geometry);
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());
122  for (auto cluster = detset.begin();
123  cluster != detset.end(); cluster++) {
124  const SiPixelCluster& aCluster = *cluster;
125  bool hasBeenSplit = false;
126  bool shouldBeSplit = false;
127  GlobalPoint cPos = det->surface().toGlobal(
128  pp->localParametersV(aCluster,
129  (*geometry->idToDetUnit(detIt->id())))[0].first);
130  GlobalPoint ppv(pv.position().x(), pv.position().y(), pv.position().z());
131  GlobalVector clusterDir = cPos - ppv;
132  for (unsigned int ji = 0; ji < cores->size(); ji++) {
133  if ((*cores)[ji].pt() > ptMin_) {
134  const reco::Candidate& jet = (*cores)[ji];
135  GlobalVector jetDir(jet.px(), jet.py(), jet.pz());
136  if (Geom::deltaR(jetDir, clusterDir) < deltaR_) {
137  // check if the cluster has to be splitted
138 
139  bool isEndCap =
140  (std::abs(cPos.z()) > 30.f); // FIXME: check detID instead!
141  float jetZOverRho = jet.momentum().Z() / jet.momentum().Rho();
142  if (isEndCap)
143  jetZOverRho = jet.momentum().Rho() / jet.momentum().Z();
144  float expSizeY =
145  std::sqrt((1.3f*1.3f) + (1.9f*1.9f) * jetZOverRho*jetZOverRho);
146  if (expSizeY < 1.f) expSizeY = 1.f;
147  float expSizeX = 1.5f;
148  if (isEndCap) {
149  expSizeX = expSizeY;
150  expSizeY = 1.5f;
151  } // in endcap col/rows are switched
152  float expCharge =
153  std::sqrt(1.08f + jetZOverRho * jetZOverRho) * centralMIPCharge_;
154 
155  if (aCluster.charge() > expCharge * chargeFracMin_ &&
156  (aCluster.sizeX() > expSizeX + 1 ||
157  aCluster.sizeY() > expSizeY + 1)) {
158  shouldBeSplit = true;
159  if (verbose)
160  std::cout << "Trying to split: charge and deltaR "
161  << aCluster.charge() << " "
162  << Geom::deltaR(jetDir, clusterDir) << " size x y "
163  << aCluster.sizeX() << " " << aCluster.sizeY()
164  << " exp. size (x,y) "
165  << expSizeX << " " << expSizeY
166  << " detid " << detIt->id() << std::endl;
167  if (verbose)
168  std::cout << "jetZOverRho=" << jetZOverRho << std::endl;
169 
170  if (split(aCluster, filler, expCharge, expSizeY, expSizeX,
171  jetZOverRho)) {
172  hasBeenSplit = true;
173  }
174  }
175  }
176  }
177  }
178  if (!hasBeenSplit) {
179  SiPixelCluster c = aCluster;
180  if (shouldBeSplit) {
181  // blowup the error if we failed to split a splittable cluster (does
182  // it ever happen)
183  c.setSplitClusterErrorX(c.sizeX() * (100.f/3.f)); // this is not really blowing up .. TODO: tune
184  c.setSplitClusterErrorY(c.sizeY() * (150.f/3.f));
185  }
186  filler.push_back(c);
187  std::push_heap(filler.begin(),filler.end(),
188  [](SiPixelCluster const & cl1,SiPixelCluster const & cl2) { return cl1.minPixelRow() < cl2.minPixelRow();});
189 
190  }
191  }// loop over clusters
192  std::sort_heap(filler.begin(),filler.end(),
193  [](SiPixelCluster const & cl1,SiPixelCluster const & cl2) { return cl1.minPixelRow() < cl2.minPixelRow();});
194 
195  } // loop over det
196  iEvent.put(std::move(output));
197 }
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:106
boost::transform_iterator< IterHelp, const_IdIter > const_iterator
const_iterator end(bool update=false) const
virtual double pz() const =0
z coordinate of momentum vector
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:122
void setSplitClusterErrorY(float erry)
edm::EDGetTokenT< edmNew::DetSetVector< SiPixelCluster > > pixelClusters_
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:460
bool split(const SiPixelCluster &aCluster, edmNew::DetSetVector< SiPixelCluster >::FastFiller &filler, float expectedADC, int sizeY, int sizeX, float jetZOverRho)
int charge() const
const GeomDetUnit * idToDetUnit(DetId) const override
Return the pointer to the GeomDetUnit corresponding to a given DetId.
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
const Point & position() const
position
Definition: Vertex.h:109
edm::EDGetTokenT< edm::View< reco::Candidate > > cores_
virtual double py() const =0
y coordinate of momentum vector
int minPixelRow() const
T sqrt(T t)
Definition: SSEVec.h:18
T z() const
Definition: PV3DBase.h:64
def pv(vc)
Definition: MetAnalyzer.py:6
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double f[11][100]
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
virtual VLocalValues localParametersV(const SiPixelCluster &cluster, const GeomDetUnit &gd) const
void setSplitClusterErrorX(float errx)
const T & get() const
Definition: EventSetup.h:55
virtual Vector momentum() const =0
spatial momentum vector
int sizeY() const
ESHandle< TrackerGeometry > geometry
Pixel cluster – collection of neighboring pixels above threshold.
HLT enums.
const GeomDet * idToDet(DetId) const override
virtual double px() const =0
x coordinate of momentum vector
int sizeX() const
T const * product() const
Definition: ESHandle.h:86
def move(src, dest)
Definition: eostools.py:510
const_iterator begin(bool update=false) const
std::multimap< float, int > JetCoreClusterSplitter::secondDistDiffScore ( const std::vector< std::vector< float > > &  distanceMap)
private

Definition at line 237 of file JetCoreClusterSplitter.cc.

References closestClusters(), and edmIntegrityCheck::d.

238  {
239  std::multimap<float, int> scores;
240  for (unsigned int j = 0; j < distanceMap.size(); j++) {
241  std::pair<float, float> d = closestClusters(distanceMap[j]);
242  scores.insert(std::pair<float, int>(d.second - d.first, j));
243  }
244  return scores;
245 }
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 247 of file JetCoreClusterSplitter.cc.

References closestClusters(), and edmIntegrityCheck::d.

Referenced by fittingSplit().

248  {
249  std::multimap<float, int> scores;
250  for (unsigned int j = 0; j < distanceMap.size(); j++) {
251  std::pair<float, float> d = closestClusters(distanceMap[j]);
252  scores.insert(std::pair<float, int>(-d.second, j));
253  }
254  return scores;
255 }
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 199 of file JetCoreClusterSplitter.cc.

References edmNew::DetSetVector< T >::FastFiller::begin(), SiPixelCluster::charge(), GetRecoTauVFromDQM_MC_cff::cl2, edmNew::DetSetVector< T >::FastFiller::end(), fittingSplit(), mps_fire::i, SiPixelCluster::minPixelRow(), and edmNew::DetSetVector< T >::FastFiller::push_back().

Referenced by produce().

202  {
203  // This function should test several configuration of splitting, and then
204  // return the one with best chi2
205 
206  std::vector<SiPixelCluster> sp =
207  fittingSplit(aCluster, expectedADC, sizeY, sizeX, jetZOverRho,
208  std::floor(aCluster.charge() / expectedADC + 0.5f));
209 
210  // for the config with best chi2
211  for (unsigned int i = 0; i < sp.size(); i++) {
212  filler.push_back(sp[i]);
213  std::push_heap(filler.begin(),filler.end(),
214  [](SiPixelCluster const & cl1,SiPixelCluster const & cl2) { return cl1.minPixelRow() < cl2.minPixelRow();});
215  }
216 
217  return (sp.size() > 0);
218 }
void push_back(data_type const &d)
int charge() const
int minPixelRow() const
std::vector< SiPixelCluster > fittingSplit(const SiPixelCluster &aCluster, float expectedADC, int sizeY, int sizeX, float jetZOverRho, unsigned int nSplitted)
Pixel cluster – collection of neighboring pixels above threshold.

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().