FastJet  3.1.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
JetConeFinderTool.cc
1 //----------------------------------------------------------------------
2 // This file distributed with FastJet has been obtained from SpartyJet
3 // v2.20.0 by Pierre-Antoine Delsart, Kurtis L. Geerlings, Joey
4 // Huston, Brian T. Martin and Chris Vermilion
5 // For details, see http://www.pa.msu.edu/~huston/SpartyJet/
6 // http://projects.hepforge.org/spartyjet/
7 //
8 // Changes from the original file are listed below.
9 //----------------------------------------------------------------------
10 
11 //*******************************************************************************
12 // Filename : JetConeFinderTool.cc
13 // Author : Ambreesh Gupta
14 // Created : Nov, 2000
15 //
16 // Jan 2004: Use CLHEP units. Use phi = (-pi,pi].
17 //*******************************************************************************
18 
19 // History of changes from the original JetConeFinder.cc file in
20 // SpartyJet v2.20
21 //
22 // 2009-01-15 Gregory Soyez <soyez@fastjet.fr>
23 //
24 // * put the code in the fastjet::atlas namespace
25 //
26 // 2009-02-14 Gregory Soyez <soyez@fastjet.fr>
27 //
28 // * imported into FastJet
29 // * removed the string name in the ctor
30 // * removed the message logs
31 // * replaced StatusCode by int
32 // * cleaned the comments
33 
34 #include <iostream>
35 
36 #include "JetConeFinderTool.hh"
37 #include "Jet.hh"
38 #include "JetDistances.hh"
39 #include "CommonUtils.hh"
40 
41 #include <vector>
42 #include <math.h>
43 
44 #include <fastjet/internal/base.hh>
45 
46 FASTJET_BEGIN_NAMESPACE
47 
48 namespace atlas {
49 
50 // set the default energy scale
51  double GeV = 1.0; //1000;
52 
53 JetConeFinderTool::JetConeFinderTool() :m_coneR(0.7)
54  , m_ptcut(0.0*GeV)
55  , m_eps(0.05)
56  , m_seedPt(2.0*GeV)
57  , m_etaMax(5.0)
58 {}
59 
60 JetConeFinderTool::~JetConeFinderTool()
61 {}
62 
63 /////////////////////////////////////////////////////////////////////////////////
64 //Execution /
65 /////////////////////////////////////////////////////////////////////////////////
66 int JetConeFinderTool::execute(jetcollection_t & theJets)
67 {
68  sort_jet_list<JetSorter_Et>(theJets);
69 
70  m_pjetV = &theJets;
71 
72  if(theJets.size()==0) return 0;
73 
74  // Initiale ctr/dctr counter for object counting.
75  m_ctr = 0;
76  m_dctr = 0;
77 
78  //////////////////////
79  // Reconstruct Jets //
80  //////////////////////
81  this->reconstruct();
82 
83  //////////////////////////
84  // ReFill JetCollection //
85  //////////////////////////
86  clear_list(theJets);
87  jetcollection_t::iterator it = m_jetOV->begin();
88  jetcollection_t::iterator itE = m_jetOV->end();
89  for(; it!=itE; ++it){
90  theJets.push_back(*it);
91  }
92 
93 
94  delete m_jetOV;
95  return 1;
96 }
97 
98 ///////////////////////////////////////////////////////////////////////////////
99 // Reconstruction algorithm specific methods /
100 ///////////////////////////////////////////////////////////////////////////////
101 
102 void
103 JetConeFinderTool::reconstruct()
104 {
105  m_jetOV = new jetcollection_t();
106 
107  jetcollection_t::iterator tItr;
108  jetcollection_t::iterator tItr_begin = m_pjetV->begin();
109  jetcollection_t::iterator tItr_end = m_pjetV->end();
110 
111  // order towers in pt
112 
113  for ( tItr=tItr_begin; tItr!=tItr_end; ++tItr ) {
114 
115  // Seed Cut
116  double tEt = (*tItr)->et();
117  if ( tEt < m_seedPt ) break;
118 
119  // Tower eta, phi
120  double etaT = (*tItr)->eta();
121  double phiT = (*tItr)->phi();
122 
123  // Iteration logic
124  bool stable = false;
125  bool inGeom = true;
126 
127  Jet* preJet;
128 
129  int count = 1;
130  do { // Iteration Loop
131 
132  // Make cone
133  preJet = calc_cone(etaT,phiT);
134  double etaC = preJet->eta();
135  double phiC = preJet->phi();
136 
137  double deta = fabs(etaT - etaC);
138  double dphi = fabs(JetDistances::deltaPhi(phiT,phiC));
139 
140  // Is Stable ?
141  if ( deta < m_eps && dphi < m_eps )
142  stable = true;
143 
144  // In Geometry ?
145  if ( fabs(etaC) > m_etaMax )
146  inGeom = false;
147 
148  etaT = etaC;
149  phiT = phiC;
150 
151  if ( !stable && inGeom ) {
152  delete preJet;
153  m_dctr +=1;
154  }
155  ++count;
156 
157  }while ( !stable && inGeom && count < 10 );
158 
159  if ( count > 9 && (!stable && inGeom) ) continue; // FIXME 9 ?
160 
161  // If iteration was succesfull -- check if this is a new jet and
162  // add it to OV.
163 
164  if ( stable && inGeom ) {
165  jetcollection_t::iterator pItr = m_jetOV->begin();
166  jetcollection_t::iterator pItrE = m_jetOV->end();
167 
168  bool newJet = true;
169  double etaT = preJet->eta();
170  double phiT = preJet->phi();
171 
172  for ( ; pItr != pItrE ; ++pItr ) {
173  double etaC = (*pItr)->eta();
174  double phiC = (*pItr)->phi();
175 
176  double deta = fabs(etaT - etaC);
177  double dphi = fabs(JetDistances::deltaPhi(phiT,phiC));
178 
179  if ( deta < 0.05 && dphi < 0.05 ) {
180  // Debugging done by Gregory Soyez:
181  //
182  // Becase of the cut on the Et difference imposed on the
183  // ordering (line 80 of Jet.hh), the ordering of the input
184  // particles in Et is not robust agains different
185  // implementations of sort (which has an undefined behaviour
186  // if 2 particles are equal). A consequence of this is that
187  // stable cone search will consider these 2 seeds in an
188  // undefined order. If the 2 resulting stable cones are too
189  // close (deta<0.05, dphi<0.05) one will be accepted and the
190  // other rejected. Which one depends on the ordering and is
191  // thus undefined. If the 2 stable cones do not have the
192  // same number of constituents this could affect the result
193  // of the clustering.
194  //
195  // The line below helps debugging these cases by printing
196  // the rejected stable cones
197  //std::cout << "rejecting " << etaT << " " << phiT << " " << preJet->et() << (*tItr)->eta() << " " << (*tItr)->phi() << " " << (*tItr)->et() << std::endl;
198  newJet = false;
199  break;
200  }
201  }
202  if ( newJet ) {
203  m_jetOV->push_back( preJet );
204  // Debugging done by Gregory Soyez:
205  //
206  // Becase of the cut on the Et difference imposed on the
207  // ordering (line 80 of Jet.hh), the ordering of the input
208  // particles in Et is not robust agains different
209  // implementations of sort (which has an undefined behaviour
210  // if 2 particles are equal). A consequence of this is that
211  // stable cone search will consider these 2 seeds in an
212  // undefined order. If the 2 resulting stable cones are too
213  // close (deta<0.05, dphi<0.05) one will be accepted and the
214  // other rejected. Which one depends on the ordering and is
215  // thus undefined. If the 2 stable cones do not have the
216  // same number of constituents this could affect the result
217  // of the clustering.
218  //
219  // The line below helps debugging these cases by printing
220  // the accepted stable cones
221  //std::cout << "accepting " << etaT << " " << phiT << " " << preJet->et() << (*tItr)->eta() << " " << (*tItr)->phi() << " " << (*tItr)->et() << std::endl;
222  }
223  else {
224  delete preJet;
225  m_dctr +=1;
226  }
227  }
228  else {
229  delete preJet;
230  m_dctr +=1;
231  }
232  }
233 }
234 
235 Jet* JetConeFinderTool::calc_cone(double eta, double phi)
236 {
237  // Create a new Jet
238  Jet* j = new Jet();
239  m_ctr +=1;
240 
241  // Add all ProtoJet within m_coneR to this Jet
242  jetcollection_t::iterator itr = m_pjetV->begin();
243  jetcollection_t::iterator itrE = m_pjetV->end();
244 
245  for ( ; itr!=itrE; ++itr ) {
246  double dR = JetDistances::deltaR(eta,phi,(*itr)->eta(),(*itr)->phi());
247  if ( dR < m_coneR ) {
248  j->addJet( (*itr) );
249  }
250  }
251 
252  return j;
253 }
254 
255 
256 
257 } // namespace atlas
258 
259 FASTJET_END_NAMESPACE