CMS 3D CMS Logo

PFMultiDepthClusterizer.cc
Go to the documentation of this file.
7 #include "Math/GenVector/VectorUtil.h"
9 
10 #include "vdt/vdtMath.h"
11 
12 #include <iterator>
13 #include <unordered_map>
14 
17 
18 public:
20 
21  ~PFMultiDepthClusterizer() override = default;
22  PFMultiDepthClusterizer(const B2DGPF&) = delete;
23  B2DGPF& operator=(const B2DGPF&) = delete;
24 
25  void update(const edm::EventSetup& es) override { _allCellsPosCalc->update(es); }
26 
28  const std::vector<bool>&,
29  reco::PFClusterCollection& outclus) override;
30 
31 private:
32  std::unique_ptr<PFCPositionCalculatorBase> _allCellsPosCalc;
33  double nSigmaEta_;
34  double nSigmaPhi_;
35 
36  class ClusterLink {
37  public:
38  ClusterLink(unsigned int i, unsigned int j, double DR, int DZ, double energy) {
39  from_ = i;
40  to_ = j;
41  linkDR_ = DR;
42  linkDZ_ = DZ;
43  linkE_ = energy;
44  }
45 
46  ~ClusterLink() = default;
47 
48  unsigned int from() const { return from_; }
49  unsigned int to() const { return to_; }
50  double dR() const { return linkDR_; }
51  int dZ() const { return linkDZ_; }
52  double energy() const { return linkE_; }
53 
54  private:
55  unsigned int from_;
56  unsigned int to_;
57  double linkDR_;
58  int linkDZ_;
59  double linkE_;
60  };
61 
62  void calculateShowerShapes(const reco::PFClusterCollection&, std::vector<double>&, std::vector<double>&);
63  std::vector<ClusterLink> link(const reco::PFClusterCollection&,
64  const std::vector<double>&,
65  const std::vector<double>&);
66  std::vector<ClusterLink> prune(std::vector<ClusterLink>&, std::vector<bool>& linkedClusters);
67 
69  unsigned int point,
70  std::vector<bool>& mask,
72  const std::vector<ClusterLink>& links);
73 
75 };
76 
78 
80  : PFClusterBuilderBase(conf, cc) {
81  if (conf.exists("allCellsPositionCalc")) {
82  const edm::ParameterSet& acConf = conf.getParameterSet("allCellsPositionCalc");
83  const std::string& algoac = acConf.getParameter<std::string>("algoName");
84  _allCellsPosCalc = PFCPositionCalculatorFactory::get()->create(algoac, acConf, cc);
85  }
86 
87  nSigmaEta_ = pow(conf.getParameter<double>("nSigmaEta"), 2);
88  nSigmaPhi_ = pow(conf.getParameter<double>("nSigmaPhi"), 2);
89 }
90 
92  const std::vector<bool>& seedable,
94  std::vector<double> etaRMS2(input.size(), 0.0);
95  std::vector<double> phiRMS2(input.size(), 0.0);
96 
97  //We need to sort the clusters for smaller to larger depth
98  // for (unsigned int i=0;i<input.size();++i)
99  // printf(" cluster%f %f \n",input[i].depth(),input[i].energy());
100 
101  //calculate cluster shapes
102  calculateShowerShapes(input, etaRMS2, phiRMS2);
103 
104  //link
105  auto&& links = link(input, etaRMS2, phiRMS2);
106  // for (const auto& link: links)
107  // printf("link %d %d %f %f\n",link.from(),link.to(),link.dR(),link.dZ());
108 
109  std::vector<bool> mask(input.size(), false);
110  std::vector<bool> linked(input.size(), false);
111 
112  //prune
113  auto&& prunedLinks = prune(links, linked);
114 
115  //printf("Pruned links\n")
116  // for (const auto& link: prunedLinks)
117  // printf("link %d %d %f %f\n",link.from(),link.to(),link.dR(),link.dZ());
118 
119  //now we need to build clusters
120  for (unsigned int i = 0; i < input.size(); ++i) {
121  //if not masked
122  if (mask[i])
123  continue;
124  //if not linked just spit it out
125  if (!linked[i]) {
126  output.push_back(input[i]);
127  // printf("Added single cluster with energy =%f \n",input[i].energy());
128  mask[i] = true;
129  continue;
130  }
131 
132  //now business: if linked and not masked gather clusters
133  reco::PFCluster cluster = input[i];
134  mask[i] = true;
135  expandCluster(cluster, i, mask, input, prunedLinks);
136  _allCellsPosCalc->calculateAndSetPosition(cluster);
137  output.push_back(cluster);
138  // printf("Added linked cluster with energy =%f\n",cluster.energy());
139  }
140 }
141 
143  std::vector<double>& etaRMS2,
144  std::vector<double>& phiRMS2) {
145  //shower shapes. here do not use the fractions
146 
147  for (unsigned int i = 0; i < clusters.size(); ++i) {
148  const reco::PFCluster& cluster = clusters[i];
149  double etaSum = 0.0;
150  double phiSum = 0.0;
151  auto const& crep = cluster.positionREP();
152  for (const auto& frac : cluster.recHitFractions()) {
153  auto const& h = *frac.recHitRef();
154  auto const& rep = h.positionREP();
155  etaSum += (frac.fraction() * h.energy()) * std::abs(rep.eta() - crep.eta());
156  phiSum += (frac.fraction() * h.energy()) * std::abs(deltaPhi(rep.phi(), crep.phi()));
157  }
158  //protection for single line : assign ~ tower
159  etaRMS2[i] = std::max(etaSum / cluster.energy(), 0.1);
160  etaRMS2[i] *= etaRMS2[i];
161  phiRMS2[i] = std::max(phiSum / cluster.energy(), 0.1);
162  phiRMS2[i] *= phiRMS2[i];
163  }
164 }
165 
166 std::vector<PFMultiDepthClusterizer::ClusterLink> PFMultiDepthClusterizer::link(
167  const reco::PFClusterCollection& clusters, const std::vector<double>& etaRMS2, const std::vector<double>& phiRMS2) {
168  std::vector<ClusterLink> links;
169  //loop on all pairs
170  for (unsigned int i = 0; i < clusters.size(); ++i)
171  for (unsigned int j = 0; j < clusters.size(); ++j) {
172  if (i == j)
173  continue;
174 
175  const reco::PFCluster& cluster1 = clusters[i];
176  const reco::PFCluster& cluster2 = clusters[j];
177 
178  // PFCluster depth stored as double but HCAL layer clusters have integral depths only
179  auto dz = (static_cast<int>(cluster2.depth()) - static_cast<int>(cluster1.depth()));
180 
181  //Do not link at the same layer and only link inside out!
182  if (dz <= 0)
183  continue;
184 
185  auto const& crep1 = cluster1.positionREP();
186  auto const& crep2 = cluster2.positionREP();
187 
188  auto deta = crep1.eta() - crep2.eta();
189  deta = deta * deta / (etaRMS2[i] + etaRMS2[j]);
190  auto dphi = deltaPhi(crep1.phi(), crep2.phi());
191  dphi = dphi * dphi / (phiRMS2[i] + phiRMS2[j]);
192 
193  // printf("Testing Link %d -> %d (%f %f %f %f ) \n",i,j,deta,dphi,cluster1.position().Eta()-cluster2.position().Eta(),deltaPhi(cluster1.position().Phi(),cluster2.position().Phi()));
194 
195  if ((deta < nSigmaEta_) & (dphi < nSigmaPhi_))
196  links.push_back(ClusterLink(i, j, deta + dphi, std::abs(dz), cluster1.energy() + cluster2.energy()));
197  }
198 
199  return links;
200 }
201 
202 std::vector<PFMultiDepthClusterizer::ClusterLink> PFMultiDepthClusterizer::prune(std::vector<ClusterLink>& links,
203  std::vector<bool>& linkedClusters) {
204  std::vector<ClusterLink> goodLinks;
205  std::vector<bool> mask(links.size(), false);
206  if (links.empty())
207  return goodLinks;
208 
209  for (unsigned int i = 0; i < links.size() - 1; ++i) {
210  if (mask[i])
211  continue;
212  for (unsigned int j = i + 1; j < links.size(); ++j) {
213  if (mask[j])
214  continue;
215 
216  const ClusterLink& link1 = links[i];
217  const ClusterLink& link2 = links[j];
218 
219  if (link1.to() == link2.to()) { //found two links going to the same spot,kill one
220  //first prefer nearby layers
221  if (link1.dZ() < link2.dZ()) {
222  mask[j] = true;
223  } else if (link1.dZ() > link2.dZ()) {
224  mask[i] = true;
225  } else { //if same layer-pick based on transverse compatibility
226  if (link1.dR() < link2.dR()) {
227  mask[j] = true;
228  } else if (link1.dR() > link2.dR()) {
229  mask[i] = true;
230  } else {
231  //same distance as well -> can happen !!!!! Pick the highest SUME
232  if (link1.energy() < link2.energy())
233  mask[i] = true;
234  else
235  mask[j] = true;
236  }
237  }
238  }
239  }
240  }
241 
242  for (unsigned int i = 0; i < links.size(); ++i) {
243  if (mask[i])
244  continue;
245  goodLinks.push_back(links[i]);
246  linkedClusters[links[i].from()] = true;
247  linkedClusters[links[i].to()] = true;
248  }
249 
250  return goodLinks;
251 }
252 
254  double e1 = 0.0;
255  double e2 = 0.0;
256 
257  //find seeds
258  for (const auto& fraction : main.recHitFractions())
259  if (fraction.recHitRef()->detId() == main.seed()) {
260  e1 = fraction.recHitRef()->energy();
261  }
262 
263  for (const auto& fraction : added.recHitFractions()) {
264  main.addRecHitFraction(fraction);
265  if (fraction.recHitRef()->detId() == added.seed()) {
266  e2 = fraction.recHitRef()->energy();
267  }
268  }
269  if (e2 > e1)
270  main.setSeed(added.seed());
271 }
272 
274  unsigned int point,
275  std::vector<bool>& mask,
277  const std::vector<ClusterLink>& links) {
278  for (const auto& link : links) {
279  if (link.from() == point) {
280  //found link that starts from this guy if not masked absorb
281  if (!mask[link.from()]) {
282  absorbCluster(cluster, clusters[link.from()]);
283  mask[link.from()] = true;
284  }
285 
286  if (!mask[link.to()]) {
287  absorbCluster(cluster, clusters[link.to()]);
288  mask[link.to()] = true;
289  expandCluster(cluster, link.to(), mask, clusters, links);
290  }
291  }
292  if (link.to() == point) {
293  //found link that starts from this guy if not masked absorb
294  if (!mask[link.to()]) {
295  absorbCluster(cluster, clusters[link.to()]);
296  mask[link.to()] = true;
297  }
298 
299  if (!mask[link.from()]) {
300  absorbCluster(cluster, clusters[link.from()]);
301  mask[link.from()] = true;
302  expandCluster(cluster, link.from(), mask, clusters, links);
303  }
304  }
305  }
306 }
void calculateShowerShapes(const reco::PFClusterCollection &, std::vector< double > &, std::vector< double > &)
T getParameter(std::string const &) const
Definition: ParameterSet.h:307
Particle flow cluster, see clustering algorithm in PFClusterAlgo.
Definition: PFCluster.h:42
DetId seed() const
return DetId of seed
Definition: CaloCluster.h:219
bool exists(std::string const &parameterName) const
checks if a parameter exists
ParameterSet const & getParameterSet(std::string const &) const
const std::vector< reco::PFRecHitFraction > & recHitFractions() const
vector of rechit fractions
Definition: PFCluster.h:65
B2DGPF & operator=(const B2DGPF &)=delete
std::vector< ClusterLink > prune(std::vector< ClusterLink > &, std::vector< bool > &linkedClusters)
const REPPoint & positionREP() const
cluster position: rho, eta, phi
Definition: PFCluster.h:92
constexpr uint32_t mask
Definition: gpuClustering.h:24
static std::string const input
Definition: EdmProvDump.cc:47
double depth() const
cluster depth
Definition: PFCluster.h:82
std::unique_ptr< PFCPositionCalculatorBase > _allCellsPosCalc
void absorbCluster(reco::PFCluster &, const reco::PFCluster &)
double energy() const
cluster energy
Definition: PFCluster.h:74
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::vector< ClusterLink > link(const reco::PFClusterCollection &, const std::vector< double > &, const std::vector< double > &)
rep
Definition: cuy.py:1189
void update(const edm::EventSetup &es) override
PFMultiDepthClusterizer(const edm::ParameterSet &conf, edm::ConsumesCollector &cc)
~PFMultiDepthClusterizer() override=default
Definition: main.py:1
std::vector< PFCluster > PFClusterCollection
collection of PFCluster objects
Definition: PFClusterFwd.h:9
#define DEFINE_EDM_PLUGIN(factory, type, name)
#define get
The Signals That Services Can Subscribe To This is based on ActivityRegistry h
Helper function to determine trigger accepts.
Definition: Activities.doc:4
void buildClusters(const reco::PFClusterCollection &, const std::vector< bool > &, reco::PFClusterCollection &outclus) override
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5
PFMultiDepthClusterizer B2DGPF
void expandCluster(reco::PFCluster &, unsigned int point, std::vector< bool > &mask, const reco::PFClusterCollection &, const std::vector< ClusterLink > &links)