CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
List of all members | Public Types | Public Member Functions | Private Member Functions | Private Attributes | Static Private Attributes | Friends
PFHcalSuperClusterAlgo Class Reference

Algorithm for particle flow superclustering. More...

#include <PFHcalSuperClusterAlgo.h>

Public Types

typedef edm::Handle
< reco::PFClusterCollection
PFClusterHandle
 
typedef edm::Ptr< reco::PFClusterPFClusterPtr
 

Public Member Functions

std::pair< double, double > calculatePosition (const reco::PFCluster &cluster)
 recalculate eta-phi position of clusters More...
 
std::pair< double, double > calculateWidths (const reco::PFCluster &cluster)
 calculate eta-phi widths of clusters More...
 
const edm::PtrVector
< reco::PFCluster > & 
clusters () const
 
std::auto_ptr< std::vector
< reco::PFCluster > > & 
clusters ()
 setters ----------------------------------------------——— More...
 
void doClustering (const reco::PFClusterCollection &clusters, const reco::PFClusterCollection &clustersHO)
 perform clustering More...
 
void doClustering (const PFClusterHandle &clustersHandle, const PFClusterHandle &clustersHOHandle)
 perform clustering in full framework More...
 
void enableDebugging (bool debug)
 enable/disable debugging More...
 
 PFHcalSuperClusterAlgo ()
 constructor More...
 
std::auto_ptr< std::vector
< reco::PFSuperCluster > > & 
superClusters ()
 
void write ()
 
virtual ~PFHcalSuperClusterAlgo ()
 destructor More...
 

Private Member Functions

void doClusteringWorker (const reco::PFClusterCollection &clusters, const reco::PFClusterCollection &clustersHO)
 perform clustering More...
 

Private Attributes

edm::PtrVector< reco::PFClusterclusters_
 
PFClusterHandle clustersHandle_
 
PFClusterHandle clustersHOHandle_
 
bool debug_
 debugging on/off More...
 
std::auto_ptr< std::vector
< reco::PFCluster > > 
pfClusters_
 particle flow clusters More...
 
std::auto_ptr< std::vector
< reco::PFSuperCluster > > 
pfSuperClusters_
 particle flow superclusters More...
 
reco::PFSuperCluster SuperCluster_
 

Static Private Attributes

static unsigned prodNum_ = 1
 product number More...
 

Friends

std::ostream & operator<< (std::ostream &out, const PFHcalSuperClusterAlgo &algo)
 

Detailed Description

Algorithm for particle flow superclustering.

This class takes as an input a map of pointers to PFCluster's, and creates PFSuperCluster's from these clusters.

Author
Chris Tully
Date
July 2012

Definition at line 27 of file PFHcalSuperClusterAlgo.h.

Member Typedef Documentation

Definition at line 41 of file PFHcalSuperClusterAlgo.h.

Definition at line 42 of file PFHcalSuperClusterAlgo.h.

Constructor & Destructor Documentation

PFHcalSuperClusterAlgo::PFHcalSuperClusterAlgo ( )

constructor

Definition at line 26 of file PFHcalSuperClusterAlgo.cc.

26  :
27  pfClusters_( new vector<reco::PFCluster> ),
28  pfSuperClusters_( new vector<reco::PFSuperCluster> ),
29  debug_(false)
30 {
31 
32 }
bool debug_
debugging on/off
std::auto_ptr< std::vector< reco::PFCluster > > pfClusters_
particle flow clusters
std::auto_ptr< std::vector< reco::PFSuperCluster > > pfSuperClusters_
particle flow superclusters
virtual PFHcalSuperClusterAlgo::~PFHcalSuperClusterAlgo ( )
inlinevirtual

destructor

Definition at line 35 of file PFHcalSuperClusterAlgo.h.

35 {;}

Member Function Documentation

std::pair< double, double > PFHcalSuperClusterAlgo::calculatePosition ( const reco::PFCluster cluster)

recalculate eta-phi position of clusters

Definition at line 60 of file PFHcalSuperClusterAlgo.cc.

References cuy::denominator, reco::PFCluster::energy(), create_public_lumi_plots::log, max(), Geom::pi(), reco::PFCluster::recHitFractions(), Geom::twoPi(), and w().

Referenced by calculateWidths(), and doClusteringWorker().

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 }
list denominator
Definition: cuy.py:484
const T & max(const T &a, const T &b)
double energy() const
cluster energy
Definition: PFCluster.h:72
double pi()
Definition: Pi.h:31
double twoPi()
Definition: Pi.h:32
const std::vector< reco::PFRecHitFraction > & recHitFractions() const
vector of rechit fractions
Definition: PFCluster.h:62
T w() const
std::pair< double, double > PFHcalSuperClusterAlgo::calculateWidths ( const reco::PFCluster cluster)

calculate eta-phi widths of clusters

Definition at line 94 of file PFHcalSuperClusterAlgo.cc.

References funct::abs(), calculatePosition(), cuy::denominator, dPhi(), reco::PFCluster::energy(), create_public_lumi_plots::log, max(), Geom::pi(), reco::PFCluster::recHitFractions(), mathSSE::sqrt(), Geom::twoPi(), and w().

Referenced by doClusteringWorker().

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 }
list denominator
Definition: cuy.py:484
double dPhi(double phi1, double phi2)
Definition: JetUtil.h:30
const T & max(const T &a, const T &b)
T sqrt(T t)
Definition: SSEVec.h:48
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double energy() const
cluster energy
Definition: PFCluster.h:72
std::pair< double, double > calculatePosition(const reco::PFCluster &cluster)
recalculate eta-phi position of clusters
double pi()
Definition: Pi.h:31
double twoPi()
Definition: Pi.h:32
const std::vector< reco::PFRecHitFraction > & recHitFractions() const
vector of rechit fractions
Definition: PFCluster.h:62
T w() const
const edm::PtrVector<reco::PFCluster>& PFHcalSuperClusterAlgo::clusters ( ) const
inline

Definition at line 43 of file PFHcalSuperClusterAlgo.h.

References clusters_.

Referenced by doClustering().

44  {return clusters_; }
edm::PtrVector< reco::PFCluster > clusters_
std::auto_ptr< std::vector< reco::PFCluster > >& PFHcalSuperClusterAlgo::clusters ( )
inline

setters ----------------------------------------------———

getters -------------------------------------------------——

Returns
particle flow clusters

Definition at line 69 of file PFHcalSuperClusterAlgo.h.

References pfClusters_.

70  {return pfClusters_;}
std::auto_ptr< std::vector< reco::PFCluster > > pfClusters_
particle flow clusters
void PFHcalSuperClusterAlgo::doClustering ( const reco::PFClusterCollection clusters,
const reco::PFClusterCollection clustersHO 
)

perform clustering

Definition at line 49 of file PFHcalSuperClusterAlgo.cc.

References edm::HandleBase::clear(), clustersHandle_, clustersHOHandle_, and doClusteringWorker().

49  {
50  // using clusters without a Handle, clear to avoid a stale member
53 
54  // perform clustering
55  doClusteringWorker( clusters, clustersHO );
56 }
void doClusteringWorker(const reco::PFClusterCollection &clusters, const reco::PFClusterCollection &clustersHO)
perform clustering
const edm::PtrVector< reco::PFCluster > & clusters() const
void PFHcalSuperClusterAlgo::doClustering ( const PFClusterHandle clustersHandle,
const PFClusterHandle clustersHOHandle 
)

perform clustering in full framework

Definition at line 37 of file PFHcalSuperClusterAlgo.cc.

References clusters(), clustersHandle_, clustersHOHandle_, and doClusteringWorker().

37  {
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 }
void doClusteringWorker(const reco::PFClusterCollection &clusters, const reco::PFClusterCollection &clustersHO)
perform clustering
std::vector< PFCluster > PFClusterCollection
collection of PFCluster objects
Definition: PFClusterFwd.h:9
const edm::PtrVector< reco::PFCluster > & clusters() const
void PFHcalSuperClusterAlgo::doClusteringWorker ( const reco::PFClusterCollection clusters,
const reco::PFClusterCollection clustersHO 
)
private

perform clustering

Definition at line 143 of file PFHcalSuperClusterAlgo.cc.

References funct::abs(), calculatePosition(), calculateWidths(), edm::PtrVectorBase::clear(), clusters_, clustersHandle_, gather_cfg::cout, SiPixelRawToDigiRegional_cfi::deltaPhi, deltaR(), dPhi(), PFRecoTauDiscriminationAgainstElectronDeadECAL_cfi::dR, first, PFLayer::HCAL_BARREL1, PFLayer::HCAL_BARREL2, PFLayer::HCAL_ENDCAP, init, reco::PFHcalSuperClusterInit::initialize(), pfClusters_, pfSuperClusters_, edm::PtrVector< T >::push_back(), edm::PtrVectorBase::size(), mathSSE::sqrt(), w2, w3, and w4.

Referenced by doClustering().

143  {
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);
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 }
size_type size() const
Size of the RefVector.
Definition: PtrVectorBase.h:73
common ppss p3p6s2 common epss epspn46 common const1 w2
Definition: inclppp.h:1
Particle flow cluster, see clustering algorithm in PFClusterAlgo.
Definition: PFCluster.h:43
void push_back(Ptr< T > const &iPtr)
Definition: PtrVector.h:138
int init
Definition: HydjetWrapper.h:62
Particle flow cluster, see clustering algorithm in PFSuperClusterAlgo.
edm::PtrVector< reco::PFCluster > clusters_
common ppss p3p6s2 common epss epspn46 common const1 w4
Definition: inclppp.h:1
double dPhi(double phi1, double phi2)
Definition: JetUtil.h:30
edm::Ptr< reco::PFCluster > PFClusterPtr
T sqrt(T t)
Definition: SSEVec.h:48
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::pair< double, double > calculateWidths(const reco::PFCluster &cluster)
calculate eta-phi widths of clusters
bool first
Definition: L1TdeRCT.cc:79
void initialize(reco::PFSuperCluster &supercluster, edm::PtrVector< reco::PFCluster > const &clusters)
double deltaR(double eta1, double eta2, double phi1, double phi2)
Definition: TreeUtility.cc:17
std::auto_ptr< std::vector< reco::PFCluster > > pfClusters_
particle flow clusters
std::pair< double, double > calculatePosition(const reco::PFCluster &cluster)
recalculate eta-phi position of clusters
common ppss p3p6s2 common epss epspn46 common const1 w3
Definition: inclppp.h:1
void clear()
Clear the PtrVector.
Definition: PtrVectorBase.h:79
tuple cout
Definition: gather_cfg.py:121
const edm::PtrVector< reco::PFCluster > & clusters() const
std::auto_ptr< std::vector< reco::PFSuperCluster > > pfSuperClusters_
particle flow superclusters
void PFHcalSuperClusterAlgo::enableDebugging ( bool  debug)
inline

enable/disable debugging

Definition at line 38 of file PFHcalSuperClusterAlgo.h.

References debug, and debug_.

38 { debug_ = debug;}
bool debug_
debugging on/off
#define debug
Definition: HDRShower.cc:19
std::auto_ptr< std::vector< reco::PFSuperCluster > >& PFHcalSuperClusterAlgo::superClusters ( )
inline
Returns
particle flow superclusters

Definition at line 73 of file PFHcalSuperClusterAlgo.h.

References pfSuperClusters_.

74  {return pfSuperClusters_;}
std::auto_ptr< std::vector< reco::PFSuperCluster > > pfSuperClusters_
particle flow superclusters
void PFHcalSuperClusterAlgo::write ( void  )

Definition at line 35 of file PFHcalSuperClusterAlgo.cc.

35  {
36 }

Friends And Related Function Documentation

std::ostream& operator<< ( std::ostream &  out,
const PFHcalSuperClusterAlgo algo 
)
friend

Member Data Documentation

edm::PtrVector< reco::PFCluster > PFHcalSuperClusterAlgo::clusters_
private

Definition at line 93 of file PFHcalSuperClusterAlgo.h.

Referenced by clusters(), and doClusteringWorker().

PFClusterHandle PFHcalSuperClusterAlgo::clustersHandle_
private

Definition at line 91 of file PFHcalSuperClusterAlgo.h.

Referenced by doClustering(), and doClusteringWorker().

PFClusterHandle PFHcalSuperClusterAlgo::clustersHOHandle_
private

Definition at line 92 of file PFHcalSuperClusterAlgo.h.

Referenced by doClustering().

bool PFHcalSuperClusterAlgo::debug_
private

debugging on/off

Definition at line 97 of file PFHcalSuperClusterAlgo.h.

Referenced by enableDebugging().

std::auto_ptr< std::vector<reco::PFCluster> > PFHcalSuperClusterAlgo::pfClusters_
private

particle flow clusters

Definition at line 85 of file PFHcalSuperClusterAlgo.h.

Referenced by clusters(), doClusteringWorker(), and operator<<().

std::auto_ptr< std::vector<reco::PFSuperCluster> > PFHcalSuperClusterAlgo::pfSuperClusters_
private

particle flow superclusters

Definition at line 88 of file PFHcalSuperClusterAlgo.h.

Referenced by doClusteringWorker(), operator<<(), and superClusters().

unsigned PFHcalSuperClusterAlgo::prodNum_ = 1
staticprivate

product number

Definition at line 100 of file PFHcalSuperClusterAlgo.h.

reco::PFSuperCluster PFHcalSuperClusterAlgo::SuperCluster_
private

Definition at line 90 of file PFHcalSuperClusterAlgo.h.