CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
PFHcalSuperClusterAlgo.cc
Go to the documentation of this file.
1 #include "RecoParticleFlow/PFClusterProducer/interface/PFHcalSuperClusterAlgo.h"
2 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
3 #include "DataFormats/Math/interface/deltaR.h"
4 #include "DataFormats/Math/interface/deltaPhi.h"
5 #include "DataFormats/DetId/interface/DetId.h"
6 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
7 #include "DataFormats/GeometryVector/interface/Pi.h"
8 #include "DataFormats/Common/interface/PtrVector.h"
9 #include "Math/GenVector/VectorUtil.h"
10 #include "RecoParticleFlow/PFClusterProducer/interface/PFHcalSuperClusterInit.h"
11 
12 #include "FWCore/MessageLogger/interface/MessageLogger.h"
13 
14 #include <stdexcept>
15 #include <string>
16 #include <sstream>
17 
18 using namespace std;
19 using namespace reco;
20 
21 unsigned PFHcalSuperClusterAlgo::prodNum_ = 1;
22 
23 //for debug only
24 //#define PFLOW_DEBUG
25 
26 PFHcalSuperClusterAlgo::PFHcalSuperClusterAlgo() :
27  pfClusters_( new vector<reco::PFCluster> ),
28  pfSuperClusters_( new vector<reco::PFSuperCluster> ),
29  debug_(false)
30 {
31 
32 }
33 
34 void
35 PFHcalSuperClusterAlgo::write() {
36 }
37 void PFHcalSuperClusterAlgo::doClustering( const PFClusterHandle& clustersHandle, const PFClusterHandle& clustersHOHandle ) {
38  const reco::PFClusterCollection& clusters = * clustersHandle;
39  const reco::PFClusterCollection& clustersHO = * clustersHOHandle;
40 
41  // cache the Handle to the clusters
42  clustersHandle_ = clustersHandle;
43  clustersHOHandle_ = clustersHOHandle;
44 
45  // perform clustering
46  doClusteringWorker( clusters, clustersHO );
47 }
48 
49 void PFHcalSuperClusterAlgo::doClustering( const reco::PFClusterCollection& clusters, const reco::PFClusterCollection& clustersHO ) {
50  // using clusters without a Handle, clear to avoid a stale member
51  clustersHandle_.clear();
52  clustersHOHandle_.clear();
53 
54  // perform clustering
55  doClusteringWorker( clusters, clustersHO );
56 }
57 
58 
59 // calculate cluster position: Rachel Myers, July 2012
60 std::pair<double, double> PFHcalSuperClusterAlgo::calculatePosition(const reco::PFCluster& cluster)
61 {
62  double numeratorEta = 0.0;
63  double numeratorPhi = 0.0;
64  double denominator = 0.0;
65  double posEta = 0.0;
66  double posPhi = 0.0;
67  double w0_ = 4.2;
68  const double clusterEnergy = cluster.energy();
69  if (cluster.energy()>0.0 && cluster.recHitFractions().size()>0) {
70  const std::vector <reco::PFRecHitFraction >& pfhitsandfracs = cluster.recHitFractions();
71  for (std::vector<reco::PFRecHitFraction>::const_iterator it = pfhitsandfracs.begin(); it != pfhitsandfracs.end(); ++it) {
72  const reco::PFRecHitRef rechit = it->recHitRef();
73  double hitEta = rechit->positionREP().Eta();
74  double hitPhi = rechit->positionREP().Phi();
75  //Change into the -pi to +pi angular range
76  while (hitPhi > +Geom::pi()) { hitPhi -= Geom::twoPi(); }
77  while (hitPhi < -Geom::pi()) { hitPhi += Geom::twoPi(); }
78  double hitEnergy = rechit->energy();
79  const double w = std::max(0.0, w0_ + log(hitEnergy / clusterEnergy));
80  denominator += w;
81  numeratorEta += w*hitEta;
82  numeratorPhi += w*hitPhi;
83  }
84  posEta = numeratorEta/denominator;
85  posPhi = numeratorPhi/denominator;
86  }
87 
88  pair<double, double> posEtaPhi(posEta,posPhi);
89 
90  return posEtaPhi;
91 }
92 
93 // calculate cluster width: Rachel Myers, July 2012
94 std::pair<double, double> PFHcalSuperClusterAlgo::calculateWidths(const reco::PFCluster& cluster)
95 {
96  double numeratorEtaEta = 0;
97  // double numeratorEtaPhi = 0;
98  double numeratorPhiPhi = 0;
99  double denominator = 0;
100  double widthEta = 0.0;
101  double widthPhi = 0.0;
102 
103  double w0_ = 4.2;
104 
105  const double clusterEta = calculatePosition(cluster).first;
106  const double clusterPhi = calculatePosition(cluster).second;
107  const double clusterEnergy = cluster.energy();
108  if(cluster.energy()>0.0 && cluster.recHitFractions().size()>0) {
109  double hitEta, hitPhi, hitEnergy, dEta, dPhi;
110  const std::vector< reco::PFRecHitFraction >& pfhitsandfracs = cluster.recHitFractions();
111  for (std::vector<reco::PFRecHitFraction>::const_iterator it = pfhitsandfracs.begin(); it != pfhitsandfracs.end(); ++it) {
112  const reco::PFRecHitRef rechit = it->recHitRef();
113  // rechit->calculatePositionREP();
114  hitEta = rechit->positionREP().Eta();
115  hitPhi = rechit->positionREP().Phi();
116  hitEnergy = rechit->energy();
117  dEta = hitEta - clusterEta;
118  dPhi = hitPhi - clusterPhi;
119 
120  while (dPhi > +Geom::pi()) { dPhi -= Geom::twoPi(); }
121  while (dPhi < -Geom::pi()) { dPhi += Geom::twoPi(); }
122 
123 
124  const double w = std::max(0.0, w0_ + log(hitEnergy / clusterEnergy));
125 
126  denominator += w;
127  numeratorEtaEta += w * dEta * dEta;
128  // numeratorEtaPhi += w * dEta * dPhi;
129  numeratorPhiPhi += w * dPhi * dPhi;
130  }
131 
132  double covEtaEta = numeratorEtaEta / denominator;
133  // double covEtaPhi_ = numeratorEtaPhi / denominator;
134  double covPhiPhi = numeratorPhiPhi / denominator;
135 
136  widthEta = sqrt(abs(covEtaEta));
137  widthPhi = sqrt(abs(covPhiPhi));
138  }
139  pair<double, double> widthEtaPhi(widthEta,widthPhi);
140  return widthEtaPhi;
141 }
142 // do clustering with new widths, positions, parameters, merging conditions: Rachel Myers, July 2012
143 void PFHcalSuperClusterAlgo::doClusteringWorker( const reco::PFClusterCollection& clusters, const reco::PFClusterCollection& clustersHO ) {
144 
145  double dRcut=0.17;
146  double dEtacut = 0.0;
147  double dPhicut = 0.0;
148  double etaScale = 1.0;
149  double phiScale = 0.5;
150  // double dRcut=0.30;
151 
152  if ( pfClusters_.get() )
153  pfClusters_->clear();
154  else
155  pfClusters_.reset( new std::vector<reco::PFCluster> );
156 
157  if ( pfSuperClusters_.get() )
158  pfSuperClusters_->clear();
159  else
160  pfSuperClusters_.reset( new std::vector<reco::PFSuperCluster> );
161 
162  // compute cluster depth index
163  std::vector< unsigned > clusterdepth(clusters.size());
164  // cout << " clusters: " << clusters.size() <<endl;
165  int mclustersHB=0;
166  int mclustersHE=0;
167  for (unsigned short ic=0; ic<clusters.size();++ic) {
168  if( clusters[ic].layer() == PFLayer::HCAL_BARREL1) mclustersHB++;
169  if( clusters[ic].layer() == PFLayer::HCAL_ENDCAP) mclustersHE++;
170  }
171  for (unsigned short ic=0; ic<clusters.size();++ic)
172  {
173  if( clusters[ic].layer() == PFLayer::HCAL_BARREL1
174  || clusters[ic].layer() == PFLayer::HCAL_ENDCAP ) { //Hcal case
175 
176  const std::vector< std::pair<DetId, float> > & hitsandfracs =
177  clusters[ic].hitsAndFractions();
178  unsigned clusterdepthfirst=0;
179  for(unsigned ihandf=0; ihandf<hitsandfracs.size(); ihandf++) {
180  unsigned depth = ((HcalDetId)hitsandfracs[ihandf].first).depth();
181  // cout << " depth parameter from clustering: " << depth <<endl;
182  clusterdepth[ic] = depth;
183  if( ihandf>0 && depth!=clusterdepthfirst) cout << " Problem with cluster depth: " << depth << " not equal to " << clusterdepthfirst <<endl;
184  else if( ihandf==0 ) clusterdepthfirst = depth;
185  }
186  // delete hitsandfracs;
187  }
188  }
189  std::vector< unsigned > clusterdepthHO(clustersHO.size());
190  // cout << " HO clusters: " << clustersHO.size() <<endl;
191  double hcaleta1=0.0;
192  double hcalphi1=0.0;
193  double hcaleta2=0.0;
194  double hcalphi2=0.0;
195  double dR = 0.0;
196  double dEta = 0.0;
197  double dPhi = 0.0;
198  for (unsigned short ic=0; ic<clustersHO.size();++ic)
199  {
200  if( clustersHO[ic].layer() == PFLayer::HCAL_BARREL2) { //HO case
201 
202  const std::vector< std::pair<DetId, float> > & hitsandfracs =
203  clustersHO[ic].hitsAndFractions();
204  unsigned clusterdepthfirst=0;
205  for(unsigned ihandf=0; ihandf<hitsandfracs.size(); ihandf++) {
206  unsigned depth = ((HcalDetId)hitsandfracs[ihandf].first).depth();
207  // cout << " depth parameter from HO clustering: " << depth <<endl;
208  clusterdepthHO[ic] = depth;
209  if( ihandf>0 && depth!=clusterdepthfirst) cout << " Problem with HO cluster depth: " << depth << " not equal to " << clusterdepthfirst <<endl;
210  else if( ihandf==0 ) clusterdepthfirst = depth;
211  }
212  // delete hitsandfracs;
213  }
214  }
215 
216  std::vector< unsigned > imerge(clusters.size());
217  std::vector< unsigned > imergeHO(clustersHO.size());
218  std::vector< bool > lmerge(clusters.size());
219  std::vector< bool > lmergeHO(clustersHO.size());
220 
221  hcaleta1=0.0;
222  hcalphi1=0.0;
223  hcaleta2=0.0;
224  hcalphi2=0.0;
225  dR = 0.0;
226  dEta = 0.0;
227  dPhi = 0.0;
228 
229  // cout << " setting up cluster merging indices "<<endl;
230  for (unsigned short ic1=0; ic1<clusters.size();++ic1) {
231  lmerge[ic1]=false;
232  }
233  for (unsigned short ic1=0; ic1<clustersHO.size();++ic1) {
234  lmergeHO[ic1]=false;
235  }
236  for (unsigned short ic1=0; ic1<clusters.size();++ic1) {
237  if( clusterdepth[ic1]==1 ){
238  hcaleta1 = calculatePosition(clusters[ic1]).first;
239  hcalphi1 = calculatePosition(clusters[ic1]).second;
240  for (unsigned short ic2=0; ic2<clusters.size();++ic2) {
241  hcaleta2 = calculatePosition(clusters[ic2]).first;
242  hcalphi2 = calculatePosition(clusters[ic2]).second;
243  dR = deltaR( hcaleta1, hcalphi1, hcaleta2, hcalphi2 );
244  dEta = abs(hcaleta1 - hcaleta2);
245  dPhi = abs(deltaPhi(hcalphi1, hcalphi2));
246  double w1 = calculateWidths(clusters[ic1]).first;
247  if (w1 == 0) w1 = 0.087;
248  double w2 = calculateWidths(clusters[ic2]).first;
249  if (w2 == 0) w2 = 0.087;
250  double w3 = calculateWidths(clusters[ic1]).second;
251  if (w3 < 0.087) w3 = 0.087;
252  double w4 = calculateWidths(clusters[ic2]).second;
253  if (w4 < 0.087) w4 = 0.087;
254  double etawidth = sqrt(w1*w1 + w2*w2);
255  double phiwidth = sqrt(w3*w3 + w4*w4);
256  dEtacut = etaScale*etawidth;
257  dPhicut = phiScale*phiwidth;
258  if( clusterdepth[ic2]==2 ){
259  // cout << " depth 1-2 dR = " << dR <<endl;
260  if((dR<dRcut) || (dEta<dEtacut && dPhi<dPhicut)) {
261  imerge[ic2]=ic1;
262  lmerge[ic2]=true;
263  }
264  } else if( clusterdepth[ic2]==3 ){
265  // cout << " depth 1-3 dR = " << dR <<endl;
266  if((dR<dRcut) || (dEta<dEtacut && dPhi<dPhicut) ) {
267  imerge[ic2]=ic1;
268  lmerge[ic2]=true;
269  }
270  }
271  }
272  } else if( clusterdepth[ic1]==2 ){
273  hcaleta1 = calculatePosition(clusters[ic1]).first;
274  hcalphi1 = calculatePosition(clusters[ic1]).second;
275  for (unsigned short ic2=0; ic2<clusters.size();++ic2) {
276  hcaleta2 = calculatePosition(clusters[ic2]).first;
277  hcalphi2 = calculatePosition(clusters[ic2]).second;
278  dEta = abs(hcaleta1 - hcaleta2);
279  dPhi = abs(deltaPhi(hcalphi1, hcalphi2));
280  dR = deltaR( hcaleta1, hcalphi1, hcaleta2, hcalphi2 );
281  double w1 = calculateWidths(clusters[ic1]).first;
282  if (w1 == 0) w1 = 0.087;
283  double w2 = calculateWidths(clusters[ic2]).first;
284  if (w2 == 0) w2 = 0.087;
285  double w3 = calculateWidths(clusters[ic1]).second;
286  if (w3 < 0.087) w3 = 0.087;
287  double w4 = calculateWidths(clusters[ic2]).second;
288  if (w4 < 0.087) w4 = 0.087;
289  double etawidth = sqrt(w1*w1+w2*w2);
290  double phiwidth = sqrt(w3*w3+w4*w4);
291  dEtacut = etaScale*etawidth;
292  dPhicut = phiScale*phiwidth;
293  if( clusterdepth[ic2]==3 ){
294  if((dR<dRcut) || (dEta<dEtacut && dPhi<dPhicut)) {
295  imerge[ic2]=ic1;
296  lmerge[ic2]=true;
297  }
298  } else if( clusterdepth[ic2]==4 ){
299  // cout << " depth 2-4 dR = " << dR <<endl;
300  if((dR<dRcut) || (dEta<dEtacut && dPhi<dPhicut)) {
301  imerge[ic2]=ic1;
302  lmerge[ic2]=true;
303  }
304  }
305  }
306  } else if( clusterdepth[ic1]==3 ){
307  hcaleta1 = calculatePosition(clusters[ic1]).first;
308  hcalphi1 = calculatePosition(clusters[ic1]).second;
309  for (unsigned short ic2=0; ic2<clusters.size();++ic2) {
310  hcaleta2 = calculatePosition(clusters[ic2]).first;
311  hcalphi2 = calculatePosition(clusters[ic2]).second;
312  dEta = abs(hcaleta1 - hcaleta2);
313  dPhi = abs(deltaPhi(hcalphi1, hcalphi2));
314  dR = deltaR( hcaleta1, hcalphi1, hcaleta2, hcalphi2 );
315  double w1 = calculateWidths(clusters[ic1]).first;
316  if (w1 == 0) w1 = 0.087;
317  double w2 = calculateWidths(clusters[ic2]).first;
318  if (w2 == 0) w2 = 0.087;
319  double w3 = calculateWidths(clusters[ic1]).second;
320  if (w3 < 0.087) w3 = 0.087;
321  double w4 = calculateWidths(clusters[ic2]).second;
322  if (w4 < 0.087) w4 = 0.087;
323  double etawidth = sqrt(w1*w1+w2*w2);
324  double phiwidth = sqrt(w3*w3+w4*w4);
325  dEtacut = etaScale*etawidth;
326  dPhicut = phiScale*phiwidth;
327  if( clusterdepth[ic2]==4 ){
328  // cout << " depth 3-4 dR = " << dR <<endl;
329  if((dR<dRcut) || (dEta<dEtacut && dPhi<dPhicut)) {
330  imerge[ic2]=ic1;
331  lmerge[ic2]=true;
332  }
333  } else if( clusterdepth[ic2]==5 ){
334  // cout << " depth 3-5 dR = " << dR <<endl;
335  if((dR<dRcut) || (dEta<dEtacut && dPhi<dPhicut)) {
336  imerge[ic2]=ic1;
337  lmerge[ic2]=true;
338  }
339  }
340  }
341  for (unsigned short ic2=0; ic2<clustersHO.size();++ic2) {
342  hcaleta2 = calculatePosition(clustersHO[ic2]).first;
343  hcalphi2 = calculatePosition(clustersHO[ic2]).second;
344  dEta = abs(hcaleta1 - hcaleta2);
345  dPhi = abs(deltaPhi(hcalphi1, hcalphi2));
346  dR = deltaR( hcaleta1, hcalphi1, hcaleta2, hcalphi2 );
347  double w1 = calculateWidths(clusters[ic1]).first;
348  if (w1 == 0) w1 = 0.087;
349  double w2 = calculateWidths(clustersHO[ic2]).first;
350  if (w2 == 0) w2 = 0.087;
351  double w3 = calculateWidths(clusters[ic1]).second;
352  if (w3 < 0.087) w3 = 0.087;
353  double w4 = calculateWidths(clustersHO[ic2]).second;
354  if (w4 < 0.087) w4 = 0.087;
355  double etawidth = sqrt(w1*w1+w2*w2);
356  double phiwidth = sqrt(w3*w3+w4*w4);
357  dEtacut = etaScale*etawidth;
358  dPhicut = phiScale*phiwidth;
359  if( clusterdepthHO[ic2]==5 ){
360  // cout << " depth 3-HO dR = " << dR <<endl;
361  if((dR<dRcut) || (dEta<dEtacut && dPhi<dPhicut)) {
362  imergeHO[ic2]=ic1;
363  lmergeHO[ic2]=true;
364  }
365  }
366  }
367  } else if( clusterdepth[ic1]==4 ){
368  hcaleta1 = calculatePosition(clusters[ic1]).first;
369  hcalphi1 = calculatePosition(clusters[ic1]).second;
370  for (unsigned short ic2=0; ic2<clusters.size();++ic2) {
371  hcaleta2 = calculatePosition(clusters[ic2]).first;
372  hcalphi2 = calculatePosition(clusters[ic2]).second;
373  dEta = abs(hcaleta1 - hcaleta2);
374  dPhi = abs(deltaPhi(hcalphi1, hcalphi2));
375  dR = deltaR( hcaleta1, hcalphi1, hcaleta2, hcalphi2 );
376  double w1 = calculateWidths(clusters[ic1]).first;
377  if (w1 < 0.087) w1 = 0.087;
378  double w2 = calculateWidths(clusters[ic2]).first;
379  if (w2 < 0.087) w2 = 0.087;
380  double w3 = calculateWidths(clusters[ic1]).second;
381  if (w3 < 0.087) w3 = 0.087;
382  double w4 = calculateWidths(clusters[ic2]).second;
383  if (w4 < 0.087) w4 = 0.087;
384  double etawidth = sqrt(w1*w1+w2*w2);
385  double phiwidth = sqrt(w3*w3+w4*w4);
386  dEtacut = etaScale*etawidth;
387  dPhicut = phiScale*phiwidth;
388  if( clusterdepth[ic2]==5 ){
389  // cout << " depth 4-5 dR = " << dR <<endl;
390  if((dR<dRcut) || (dEta<dEtacut && dPhi<dPhicut)) {
391  imerge[ic2]=ic1;
392  lmerge[ic2]=true;
393  }
394  }
395  }
396  } else if( clusterdepth[ic1]==5 ){
397  hcaleta1 = calculatePosition(clusters[ic1]).first;
398  hcalphi1 = calculatePosition(clusters[ic1]).second;
399  }
400  }
401 
402  // start a supercluster with a depth=1 cluster, then loop on all other
403  // clusters to add to cluster list, then for each cluster to add, loop
404  // on remaining clusters to check 2nd level of addition, repeat for all layers
405 
406  // need to add HO cluster logic
407  edm::PtrVector< reco::PFCluster > mergeclusters;
408  std::vector< bool > lmergeclusters(clusters.size());
409  for (unsigned short id=0; id<4;++id) {
410  // cout << " merging with starting depth: "<<id<<endl;
411  for (unsigned short ic1=0; ic1<clusters.size();++ic1) {
412  if( clusterdepth[ic1]==(unsigned)(1+id) ){
413  if(!lmerge[ic1]) {
414  for (unsigned short ic=0; ic<clusters.size();++ic) {
415  lmergeclusters[ic]=false;
416  }
417  // mergeclusters.push_back(clusters[ic1]);
418  for (unsigned short ic2=0; ic2<clusters.size();++ic2) {
419  if( clusterdepth[ic2]==(unsigned)(2+id) ){
420  if( imerge[ic2]==ic1 && lmerge[ic2] ) {
421  // mergeclusters.push_back(clusters[ic2]);
422  lmergeclusters[ic2]=true;
423  for (unsigned short ic3=0; ic3<clusters.size();++ic3) {
424  if( clusterdepth[ic3]==(unsigned)(3+id) ){
425  if( imerge[ic3]==ic2 && lmerge[ic3] ) {
426  // mergeclusters.push_back(clusters[ic3]);
427  lmergeclusters[ic3]=true;
428  for (unsigned short ic4=0; ic4<clusters.size();++ic4) {
429  if( clusterdepth[ic4]==(unsigned)(4+id) ){
430  if( imerge[ic4]==ic3 && lmerge[ic4] ) {
431  // mergeclusters.push_back(clusters[ic4]);
432  lmergeclusters[ic4]=true;
433  for (unsigned short ic5=0; ic5<clusters.size();++ic5) {
434  if( clusterdepth[ic5]==(unsigned)(5+id) ){
435  // if( imerge[ic5]==ic4 && lmerge[ic5] ) mergeclusters.push_back(clusters[ic5]);
436  if( imerge[ic5]==ic4 && lmerge[ic5] ) lmergeclusters[ic5]=true;
437  }
438  }
439  } else if( clusterdepth[ic4]==(unsigned)(5+id) ){
440  // if( imerge[ic4]==ic3 && lmerge[ic4] ) mergeclusters.push_back(clusters[ic4]);
441  if( imerge[ic4]==ic3 && lmerge[ic4] ) lmergeclusters[ic4]=true;
442  }
443  }
444  }
445  } else if( clusterdepth[ic3]==(unsigned)(4+id) ){
446  if( imerge[ic3]==ic2 && lmerge[ic3] ) {
447  // mergeclusters.push_back(clusters[ic3]);
448  lmergeclusters[ic3]=true;
449  for (unsigned short ic4=0; ic4<clusters.size();++ic4) {
450  if( clusterdepth[ic4]==(unsigned)(5+id) ){
451  // if( imerge[ic4]==ic3 && lmerge[ic4] ) mergeclusters.push_back(clusters[ic4]);
452  if( imerge[ic4]==ic3 && lmerge[ic4] ) lmergeclusters[ic4]=true;
453  }
454  }
455  }
456  }
457  } if( clusterdepth[ic3]==(unsigned)(5+id) ){
458  // if( imerge[ic3]==ic2 && lmerge[ic3] ) mergeclusters.push_back(clusters[ic3]);
459  if( imerge[ic3]==ic2 && lmerge[ic3] ) lmergeclusters[ic3]=true;
460  }
461  }
462  }
463  } else if( clusterdepth[ic2]==(unsigned)(3+id) ){
464  if( imerge[ic2]==ic1 && lmerge[ic2] ) {
465  // mergeclusters.push_back(clusters[ic2]);
466  lmergeclusters[ic2]=true;
467  for (unsigned short ic3=0; ic3<clusters.size();++ic3) {
468  if( clusterdepth[ic3]==(unsigned)(4+id) ){
469  if( imerge[ic3]==ic2 && lmerge[ic3] ) {
470  // mergeclusters.push_back(clusters[ic3]);
471  lmergeclusters[ic3]=true;
472  for (unsigned short ic4=0; ic4<clusters.size();++ic4) {
473  if( clusterdepth[ic4]==(unsigned)(5+id) ){
474  // if( imerge[ic4]==ic3 && lmerge[ic4] ) mergeclusters.push_back(clusters[ic4]);
475  if( imerge[ic4]==ic3 && lmerge[ic4] ) lmergeclusters[ic4]=true;
476  }
477  }
478  } else if( clusterdepth[ic3]==(unsigned)(5+id) ){
479  // if( imerge[ic3]==ic2 && lmerge[ic3] ) mergeclusters.push_back(clusters[ic3]);
480  if( imerge[ic3]==ic2 && lmerge[ic3] ) lmergeclusters[ic3]=true;
481  }
482  }
483  }
484  }
485  }
486  }
487  mergeclusters.push_back( PFClusterPtr( clustersHandle_, ic1));
488  for (unsigned short ic=0; ic<clusters.size();++ic) {
489  if(lmergeclusters[ic]) {
490  mergeclusters.push_back( PFClusterPtr(clustersHandle_, ic ) );
491  }
492  }
493  if(mergeclusters.size()>0) {
494  // cout << " number of clusters to merge: " <<mergeclusters.size()<<endl;
495  reco::PFSuperCluster ipfsupercluster(mergeclusters);
496  PFHcalSuperClusterInit init;
497  init.initialize( ipfsupercluster, clusters_);
498  pfSuperClusters_->push_back(ipfsupercluster);
499  pfClusters_->push_back((reco::PFCluster)ipfsupercluster);
500  mergeclusters.clear();
501  }
502 
503  }
504  } // end of depth 1+id initiated logic
505  }
506  }
507 
508  /*
509  for (unsigned short ic=0; ic<clusters.size();++ic) {
510  mergeclusters.clear();
511  mergeclusters.push_back(clusters[ic]);
512  reco::PFSuperCluster ipfsupercluster(mergeclusters);
513  pfSuperClusters_->push_back(ipfsupercluster);
514  pfClusters_->push_back((reco::PFCluster)ipfsupercluster);
515  // pfClusters_->push_back(clusters[ic]);
516  }
517  */
518 
519  clusterdepth.clear();
520  clusterdepthHO.clear();
521  imerge.clear();
522  imergeHO.clear();
523  lmerge.clear();
524  lmergeHO.clear();
525  mergeclusters.clear();
526 
527 }
528 ostream& operator<<(ostream& out,const PFHcalSuperClusterAlgo& algo) {
529  if(!out) return out;
530  out<<"PFSuperClusterAlgo parameters : "<<endl;
531  out<<"-----------------------------------------------------"<<endl;
532 
533  out<<endl;
534  out<<algo.pfClusters_->size()<<" clusters:"<<endl;
535 
536  for(unsigned i=0; i<algo.pfClusters_->size(); i++) {
537  out<<(*algo.pfClusters_)[i]<<endl;
538 
539  if(!out) return out;
540  }
541 
542  out<<algo.pfSuperClusters_->size()<<" superclusters:"<<endl;
543 
544  for(unsigned i=0; i<algo.pfSuperClusters_->size(); i++) {
545  out<<(*algo.pfSuperClusters_)[i]<<endl;
546 
547  if(!out) return out;
548  }
549  return out;
550 }
i
int i
Definition: DBlmapReader.cc:9
edm::PtrVectorBase::size
size_type size() const
Size of the RefVector.
Definition: PtrVectorBase.h:73
w2
common ppss p3p6s2 common epss epspn46 common const1 w2
Definition: inclppp.h:1
reco::PFCluster
Particle flow cluster, see clustering algorithm in PFClusterAlgo.
Definition: PFCluster.h:43
edm::PtrVector::push_back
void push_back(Ptr< T > const &iPtr)
Definition: PtrVector.h:138
init
int init
Definition: HydjetWrapper.h:62
create_public_lumi_plots.log
log
Definition: create_public_lumi_plots.py:1102
reco::PFSuperCluster
Particle flow cluster, see clustering algorithm in PFSuperClusterAlgo.
Definition: PFSuperCluster.h:26
operator<<
std::ostream & operator<<(std::ostream &out, const ALILine &li)
Definition: ALILine.cc:187
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:105
PFHcalSuperClusterAlgo::prodNum_
static unsigned prodNum_
product number
Definition: PFHcalSuperClusterAlgo.h:100
PFHcalSuperClusterAlgo::clusters_
edm::PtrVector< reco::PFCluster > clusters_
Definition: PFHcalSuperClusterAlgo.h:93
w4
common ppss p3p6s2 common epss epspn46 common const1 w4
Definition: inclppp.h:1
cuy.denominator
list denominator
Definition: cuy.py:484
dPhi
double dPhi(double phi1, double phi2)
Definition: JetUtil.h:30
max
const T & max(const T &a, const T &b)
Definition: MaterialBudgetTrackerHistos.cc:4
PFHcalSuperClusterAlgo::PFClusterPtr
edm::Ptr< reco::PFCluster > PFClusterPtr
Definition: PFHcalSuperClusterAlgo.h:42
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
PFHcalSuperClusterAlgo::doClusteringWorker
void doClusteringWorker(const reco::PFClusterCollection &clusters, const reco::PFClusterCollection &clustersHO)
perform clustering
Definition: PFHcalSuperClusterAlgo.cc:143
PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi.dR
tuple dR
Definition: PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi.py:19
PFHcalSuperClusterAlgo::calculateWidths
std::pair< double, double > calculateWidths(const reco::PFCluster &cluster)
calculate eta-phi widths of clusters
Definition: PFHcalSuperClusterAlgo.cc:94
first
bool first
Definition: L1TdeRCT.cc:79
reco::PFCluster::energy
double energy() const
cluster energy
Definition: PFCluster.h:72
reco::PFHcalSuperClusterInit::initialize
void initialize(reco::PFSuperCluster &supercluster, edm::PtrVector< reco::PFCluster > const &clusters)
Definition: PFHcalSuperClusterInit.cc:10
dbtoconf.out
tuple out
Definition: dbtoconf.py:99
deltaR
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
PFHcalSuperClusterAlgo::pfClusters_
std::auto_ptr< std::vector< reco::PFCluster > > pfClusters_
particle flow clusters
Definition: PFHcalSuperClusterAlgo.h:85
SiPixelRawToDigiRegional_cfi.deltaPhi
tuple deltaPhi
Definition: SiPixelRawToDigiRegional_cfi.py:9
PFHcalSuperClusterAlgo::calculatePosition
std::pair< double, double > calculatePosition(const reco::PFCluster &cluster)
recalculate eta-phi position of clusters
Definition: PFHcalSuperClusterAlgo.cc:60
w3
common ppss p3p6s2 common epss epspn46 common const1 w3
Definition: inclppp.h:1
edm::PtrVectorBase::clear
void clear()
Clear the PtrVector.
Definition: PtrVectorBase.h:79
Geom::pi
double pi()
Definition: Pi.h:31
Geom::twoPi
double twoPi()
Definition: Pi.h:32
PFHcalSuperClusterAlgo
Algorithm for particle flow superclustering.
Definition: PFHcalSuperClusterAlgo.h:27
reco::PFClusterCollection
std::vector< PFCluster > PFClusterCollection
collection of PFCluster objects
Definition: PFClusterFwd.h:9
reco::PFCluster::recHitFractions
const std::vector< reco::PFRecHitFraction > & recHitFractions() const
vector of rechit fractions
Definition: PFCluster.h:62
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
w
T w() const
Definition: extBasic3DVector.h:245
edm::false
volatile std::atomic< bool > shutdown_flag false
Definition: UnixSignalHandlers.cc:22
PFHcalSuperClusterAlgo::clusters
const edm::PtrVector< reco::PFCluster > & clusters() const
Definition: PFHcalSuperClusterAlgo.h:43
PFHcalSuperClusterAlgo::doClustering
void doClustering(const reco::PFClusterCollection &clusters, const reco::PFClusterCollection &clustersHO)
perform clustering
Definition: PFHcalSuperClusterAlgo.cc:49
PFHcalSuperClusterAlgo::pfSuperClusters_
std::auto_ptr< std::vector< reco::PFSuperCluster > > pfSuperClusters_
particle flow superclusters
Definition: PFHcalSuperClusterAlgo.h:88
Generated for CMSSW Reference Manual by doxygen 1.8.5