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)
 
 ProducerBase ()
 
void registerProducts (ProducerBase *, ProductRegistry *, ModuleDescription const &)
 
std::function< void(BranchDescription
const &)> 
registrationCallback () const
 used by the fwk to register list of products More...
 
virtual ~ProducerBase ()
 
- Public Member Functions inherited from edm::EDConsumerBase
std::vector< ConsumesInfoconsumesInfo () const
 
 EDConsumerBase ()
 
ProductHolderIndexAndSkipBit indexFrom (EDGetToken, BranchType, TypeID const &) const
 
void itemsMayGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const
 
void itemsToGet (BranchType, std::vector< ProductHolderIndexAndSkipBit > &) const
 
std::vector
< ProductHolderIndexAndSkipBit >
const & 
itemsToGetFromEvent () const
 
void labelsForToken (EDGetToken iToken, Labels &oLabels) const
 
void modulesDependentUpon (const std::string &iProcessName, std::vector< const char * > &oModuleLabels) const
 
void modulesWhoseProductsAreConsumed (std::vector< ModuleDescription const * > &modules, ProductRegistry const &preg, std::map< std::string, ModuleDescription const * > const &labelsToDesc, std::string const &processName) const
 
bool registeredToConsume (ProductHolderIndex, bool, BranchType) const
 
bool registeredToConsumeMany (TypeID const &, BranchType) const
 
void updateLookup (BranchType iBranchType, ProductHolderIndexHelper const &)
 
virtual ~EDConsumerBase ()
 

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 217 of file JetCoreClusterSplitter.cc.

References i.

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

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

Definition at line 253 of file JetCoreClusterSplitter.cc.

References closestClusters(), and j.

254  {
255  std::multimap<float, int> scores;
256  for (unsigned int j = 0; j < distanceMap.size(); j++) {
257  std::pair<float, float> d = closestClusters(distanceMap[j]);
258  scores.insert(std::pair<float, int>(-d.first, j));
259  }
260  return scores;
261 }
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 263 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_steps::k, convertSQLitetoXML_cfg::output, AlCaHLTBitMon_ParallelJobs::p, SiPixelCluster::pixels(), edm::second(), secondDistScore(), findQualityFiles::size, mathSSE::sqrt(), x, and y.

Referenced by split().

265  {
266  std::vector<SiPixelCluster> output;
267 
268  unsigned int meanExp = nSplitted;
269  if (meanExp <= 1) {
270  output.push_back(aCluster);
271  return output;
272  }
273 
274  std::vector<float> clx(meanExp);
275  std::vector<float> cly(meanExp);
276  std::vector<float> cls(meanExp);
277  std::vector<float> oldclx(meanExp);
278  std::vector<float> oldcly(meanExp);
279  std::vector<SiPixelCluster::Pixel> originalpixels = aCluster.pixels();
280  std::vector<std::pair<int, SiPixelCluster::Pixel> > pixels;
281  for (unsigned int j = 0; j < originalpixels.size(); j++) {
282  int sub = originalpixels[j].adc / chargePerUnit_ * expectedADC /
284  if (sub < 1) sub = 1;
285  int perDiv = originalpixels[j].adc / sub;
286  if (verbose)
287  std::cout << "Splitting " << j << " in [ " << pixels.size() << " , "
288  << pixels.size() + sub << " ], expected numb of clusters: "
289  << meanExp << " original pixel (x,y) "
290  << originalpixels[j].x << " " << originalpixels[j].y
291  << " sub " << sub << std::endl;
292  for (int k = 0; k < sub; k++) {
293  if (k == sub - 1) perDiv = originalpixels[j].adc - perDiv * k;
294  pixels.push_back(std::make_pair(j, SiPixelCluster::Pixel(originalpixels[j].x,
295  originalpixels[j].y, perDiv)));
296  }
297  }
298  std::vector<int> clusterForPixel(pixels.size());
299  // initial values
300  for (unsigned int j = 0; j < meanExp; j++) {
301  oldclx[j] = -999;
302  oldcly[j] = -999;
303  clx[j] = originalpixels[0].x + j;
304  cly[j] = originalpixels[0].y + j;
305  cls[j] = 0;
306  }
307  bool stop = false;
308  int remainingSteps = 100;
309  while (!stop && remainingSteps > 0) {
310  remainingSteps--;
311  // Compute all distances
312  std::vector<std::vector<float> > distanceMapX(originalpixels.size(),
313  std::vector<float>(meanExp));
314  std::vector<std::vector<float> > distanceMapY(originalpixels.size(),
315  std::vector<float>(meanExp));
316  std::vector<std::vector<float> > distanceMap(originalpixels.size(),
317  std::vector<float>(meanExp));
318  for (unsigned int j = 0; j < originalpixels.size(); j++) {
319  if (verbose)
320  std::cout << "Original Pixel pos " << j << " " << pixels[j].second.x << " "
321  << pixels[j].second.y << std::endl;
322  for (unsigned int i = 0; i < meanExp; i++) {
323  distanceMapX[j][i] = 1. * originalpixels[j].x - clx[i];
324  distanceMapY[j][i] = 1. * originalpixels[j].y - cly[i];
325  float dist = 0;
326  // float sizeX=2;
327  if (std::abs(distanceMapX[j][i]) > sizeX / 2.) {
328  dist += (std::abs(distanceMapX[j][i]) - sizeX / 2. + 1) *
329  (std::abs(distanceMapX[j][i]) - sizeX / 2. + 1);
330  } else {
331  dist +=
332  distanceMapX[j][i] / sizeX * 2 * distanceMapX[j][i] / sizeX * 2;
333  }
334 
335  if (std::abs(distanceMapY[j][i]) > sizeY / 2.) {
336  dist += 1. * (std::abs(distanceMapY[j][i]) - sizeY / 2. + 1.) *
337  (std::abs(distanceMapY[j][i]) - sizeY / 2. + 1.);
338  } else {
339  dist += 1. * distanceMapY[j][i] / sizeY * 2. * distanceMapY[j][i] /
340  sizeY * 2.;
341  }
342  distanceMap[j][i] = sqrt(dist);
343  if (verbose)
344  std::cout << "Cluster " << i << " Original Pixel " << j
345  << " distances: " << distanceMapX[j][i] << " "
346  << distanceMapY[j][i] << " " << distanceMap[j][i]
347  << std::endl;
348  }
349  }
350  // Compute scores for sequential addition. The first index is the
351  // distance, in whatever metrics we use, while the second is the
352  // pixel index w.r.t which the distance is computed.
353  std::multimap<float, int> scores;
354 
355  // Using different rankings to improve convergence (as Giulio proposed)
356  scores = secondDistScore(distanceMap);
357 
358  // Iterate starting from the ones with furthest second best clusters, i.e.
359  // easy choices
360  std::vector<float> weightOfPixel(pixels.size());
361  for (std::multimap<float, int>::iterator it = scores.begin();
362  it != scores.end(); it++) {
363  int pixel_index = it->second;
364  if (verbose)
365  std::cout << "Original Pixel " << pixel_index << " with score " << it->first << std::endl;
366  // find cluster that is both close and has some charge still to assign
367  int subpixel_counter = 0;
368  for (auto subpixel = pixels.begin(); subpixel != pixels.end();
369  ++subpixel, ++subpixel_counter) {
370  if (subpixel->first > pixel_index) {
371  break;
372  } else if (subpixel->first != pixel_index) {
373  continue;
374  } else {
375  float maxEst = 0;
376  int cl = -1;
377  for (unsigned int subcluster_index = 0;
378  subcluster_index < meanExp; subcluster_index++) {
379  float nsig =
380  (cls[subcluster_index] - expectedADC) /
381  (expectedADC *
382  fractionalWidth_); // 20% uncertainty? realistic from Landau?
383  float clQest = 1. / (1. + exp(nsig)) + 1e-6; // 1./(1.+exp(x*x-3*3))
384  float clDest = 1. / (distanceMap[pixel_index][subcluster_index] + 0.05);
385 
386  if (verbose)
387  std::cout << " Q: " << clQest << " D: " << clDest << " "
388  << distanceMap[pixel_index][subcluster_index] << std::endl;
389  float est = clQest * clDest;
390  if (est > maxEst) {
391  cl = subcluster_index;
392  maxEst = est;
393  }
394  }
395  cls[cl] += subpixel->second.adc;
396  clusterForPixel[subpixel_counter] = cl;
397  weightOfPixel[subpixel_counter] = maxEst;
398  if (verbose)
399  std::cout << "Pixel weight j cl " << weightOfPixel[subpixel_counter]
400  << " " << subpixel_counter
401  << " " << cl << std::endl;
402  }
403  }
404  }
405  // Recompute cluster centers
406  stop = true;
407  for (unsigned int subcluster_index = 0;
408  subcluster_index < meanExp; subcluster_index++) {
409  if (std::abs(clx[subcluster_index] - oldclx[subcluster_index]) > 0.01)
410  stop = false; // still moving
411  if (std::abs(cly[subcluster_index] - oldcly[subcluster_index]) > 0.01)
412  stop = false;
413  oldclx[subcluster_index] = clx[subcluster_index];
414  oldcly[subcluster_index] = cly[subcluster_index];
415  clx[subcluster_index] = 0;
416  cly[subcluster_index] = 0;
417  cls[subcluster_index] = 1e-99;
418  }
419  for (unsigned int pixel_index = 0;
420  pixel_index < pixels.size(); pixel_index++) {
421  if (clusterForPixel[pixel_index] < 0) continue;
422  if (verbose)
423  std::cout << "j " << pixel_index << " x " << pixels[pixel_index].second.x << " * y "
424  << pixels[pixel_index].second.y << " * ADC "
425  << pixels[pixel_index].second.adc << " * W "
426  << weightOfPixel[pixel_index] << std::endl;
427  clx[clusterForPixel[pixel_index]] += pixels[pixel_index].second.x * pixels[pixel_index].second.adc;
428  cly[clusterForPixel[pixel_index]] += pixels[pixel_index].second.y * pixels[pixel_index].second.adc;
429  cls[clusterForPixel[pixel_index]] += pixels[pixel_index].second.adc;
430  }
431  for (unsigned int subcluster_index = 0;
432  subcluster_index < meanExp; subcluster_index++) {
433  if (cls[subcluster_index] != 0) {
434  clx[subcluster_index] /= cls[subcluster_index];
435  cly[subcluster_index] /= cls[subcluster_index];
436  }
437  if (verbose)
438  std::cout << "Center for cluster " << subcluster_index << " x,y "
439  << clx[subcluster_index] << " "
440  << cly[subcluster_index] << std::endl;
441  cls[subcluster_index] = 0;
442  }
443  }
444  if (verbose) std::cout << "maxstep " << remainingSteps << std::endl;
445  // accumulate pixel with same cl
446  std::vector<std::vector<SiPixelCluster::Pixel> > pixelsForCl(meanExp);
447  for (int cl = 0; cl < (int)meanExp; cl++) {
448  for (unsigned int j = 0; j < pixels.size(); j++) {
449  if (clusterForPixel[j] == cl and
450  pixels[j].second.adc != 0) { // for each pixel of cluster
451  // cl find the other pixels
452  // with same x,y and
453  // accumulate+reset their adc
454  for (unsigned int k = j + 1; k < pixels.size(); k++) {
455  if (pixels[k].second.adc != 0
456  and pixels[k].second.x == pixels[j].second.x
457  and pixels[k].second.y == pixels[j].second.y
458  and clusterForPixel[k] == cl) {
459  if (verbose)
460  std::cout << "Resetting all sub-pixel for location "
461  << pixels[k].second.x << ", " << pixels[k].second.y
462  << " at index " << k << " associated to cl "
463  << clusterForPixel[k] << std::endl;
464  pixels[j].second.adc += pixels[k].second.adc;
465  pixels[k].second.adc = 0;
466  }
467  }
468  for (unsigned int p = 0; p < pixels.size(); ++p)
469  if (verbose)
470  std::cout << "index, x, y, ADC: " << p << ", "
471  << pixels[p].second.x << ", " << pixels[p].second.y
472  << ", " << pixels[p].second.adc
473  << " associated to cl " << clusterForPixel[p] << std::endl
474  << "Adding pixel " << pixels[j].second.x << ", " << pixels[j].second.y
475  << " to cluster " << cl << std::endl;
476  pixelsForCl[cl].push_back(pixels[j].second);
477  }
478  }
479  }
480 
481  // std::vector<std::vector<std::vector<SiPixelCluster::PixelPos *> > >
482  //pixelMap(meanExp,std::vector<std::vector<SiPixelCluster::PixelPos *>
483  //>(512,std::vector<SiPixelCluster::Pixel *>(512,0)));
484 
485  for (int cl = 0; cl < (int)meanExp; cl++) {
486  if (verbose) std::cout << "Pixels of cl " << cl << " ";
487  for (unsigned int j = 0; j < pixelsForCl[cl].size(); j++) {
488  SiPixelCluster::PixelPos newpix(pixelsForCl[cl][j].x,
489  pixelsForCl[cl][j].y);
490  if (verbose)
491  std::cout << pixelsForCl[cl][j].x << "," << pixelsForCl[cl][j].y << "|";
492  if (j==0) {
493  output.emplace_back(newpix, pixelsForCl[cl][j].adc);
494  } else {
495  output.back().add(newpix, pixelsForCl[cl][j].adc);
496  }
497  }
498  if (verbose) std::cout << std::endl;
499  if (pixelsForCl[cl].size() > 0) {
500  if (forceXError_ > 0) output.back().setSplitClusterErrorX(forceXError_);
501  if (forceYError_ > 0) output.back().setSplitClusterErrorY(forceYError_);
502  }
503  }
504  // if(verbose) std::cout << "Weights" << std::endl;
505  // if(verbose) print(theWeights,aCluster,1);
506  // if(verbose) std::cout << "Unused charge" << std::endl;
507  // if(verbose) print(theBufferResidual,aCluster);
508 
509  return output;
510 }
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:48
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int j
Definition: DBlmapReader.cc:9
tuple cout
Definition: gather_cfg.py:121
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_, 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(), 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_25ns14e33_v1_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  std::auto_ptr<edmNew::DetSetVector<SiPixelCluster> > output(
115 
117  inputPixelClusters->begin();
118  for (; detIt != inputPixelClusters->end(); detIt++) {
120  detIt->id());
121  const edmNew::DetSet<SiPixelCluster>& detset = *detIt;
122  const GeomDet* det = geometry->idToDet(detset.id());
124  detset.begin();
125  cluster != detset.end(); cluster++) {
126  const SiPixelCluster& aCluster = *cluster;
127  bool hasBeenSplit = false;
128  bool shouldBeSplit = false;
129  GlobalPoint cPos = det->surface().toGlobal(
130  pp->localParametersV(aCluster,
131  (*geometry->idToDetUnit(detIt->id())))[0].first);
132  GlobalPoint ppv(pv.position().x(), pv.position().y(), pv.position().z());
133  GlobalVector clusterDir = cPos - ppv;
134  for (unsigned int ji = 0; ji < cores->size(); ji++) {
135  if ((*cores)[ji].pt() > ptMin_) {
136  const reco::Candidate& jet = (*cores)[ji];
137  GlobalVector jetDir(jet.px(), jet.py(), jet.pz());
138  if (Geom::deltaR(jetDir, clusterDir) < deltaR_) {
139  // check if the cluster has to be splitted
140 
141  bool isEndCap =
142  (fabs(cPos.z()) > 30); // FIXME: check detID instead!
143  float jetZOverRho = jet.momentum().Z() / jet.momentum().Rho();
144  if (isEndCap)
145  jetZOverRho = jet.momentum().Rho() / jet.momentum().Z();
146  float expSizeY =
147  fabs(sqrt(1.3 * 1.3 + 1.9 * 1.9 * jetZOverRho * jetZOverRho));
148  if (expSizeY < 1) expSizeY = 1.;
149  float expSizeX = 1.5;
150  if (isEndCap) {
151  expSizeX = expSizeY;
152  expSizeY = 1.5;
153  } // in endcap col/rows are switched
154  float expCharge =
155  sqrt(1.08 + jetZOverRho * jetZOverRho) * centralMIPCharge_;
156 
157  if (aCluster.charge() > expCharge * chargeFracMin_ &&
158  (aCluster.sizeX() > expSizeX + 1 ||
159  aCluster.sizeY() > expSizeY + 1)) {
160  shouldBeSplit = true;
161  if (verbose)
162  std::cout << "Trying to split: charge and deltaR "
163  << aCluster.charge() << " "
164  << Geom::deltaR(jetDir, clusterDir) << " size x y "
165  << aCluster.sizeX() << " " << aCluster.sizeY()
166  << " exp. size (x,y) "
167  << expSizeX << " " << expSizeY
168  << " detid " << detIt->id() << std::endl;
169  if (verbose)
170  std::cout << "jetZOverRho=" << jetZOverRho << std::endl;
171 
172  if (split(aCluster, filler, expCharge, expSizeY, expSizeX,
173  jetZOverRho)) {
174  hasBeenSplit = true;
175  }
176  }
177  }
178  }
179  }
180  if (!hasBeenSplit) {
181  SiPixelCluster c = aCluster;
182  if (shouldBeSplit) {
183  // blowup the error if we failed to split a splittable cluster (does
184  // it ever happen)
186  c.sizeX() * 100. /
187  3.); // this is not really blowing up .. TODO: tune
188  c.setSplitClusterErrorY(c.sizeY() * 150. / 3.);
189  }
190  filler.push_back(c);
191  }
192  }
193  }
194  iEvent.put(output);
195 }
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:114
float charge() const
boost::transform_iterator< IterHelp, const_IdIter > const_iterator
const_iterator end(bool update=false) const
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:40
edm::EDGetTokenT< edm::View< reco::Candidate > > cores_
virtual Vector momentum() const =0
spatial momentum vector
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:115
T sqrt(T t)
Definition: SSEVec.h:48
T z() const
Definition: PV3DBase.h:64
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:55
int sizeY() const
ESHandle< TrackerGeometry > geometry
Pixel cluster – collection of neighboring pixels above threshold.
tuple cout
Definition: gather_cfg.py:121
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 233 of file JetCoreClusterSplitter.cc.

References closestClusters(), and j.

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

References closestClusters(), and j.

Referenced by fittingSplit().

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

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

Referenced by produce().

200  {
201  // This function should test several configuration of splitting, and then
202  // return the one with best chi2
203 
204  std::vector<SiPixelCluster> sp =
205  fittingSplit(aCluster, expectedADC, sizeY, sizeX, jetZOverRho,
206  floor(aCluster.charge() / expectedADC + 0.5));
207 
208  // for the config with best chi2
209  for (unsigned int i = 0; i < sp.size(); i++) {
210  filler.push_back(sp[i]);
211  }
212 
213  return (sp.size() > 0);
214 }
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)

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