Fork me on GitHub

source: git/external/fastjet/plugins/ATLASCone/JetConeFinderTool.cc@ 7b0e00c

Last change on this file since 7b0e00c was d7d2da3, checked in by pavel <pavel@…>, 12 years ago

move branches/ModularDelphes to trunk

  • Property mode set to 100644
File size: 7.3 KB
RevLine 
[d7d2da3]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
46FASTJET_BEGIN_NAMESPACE
47
48namespace atlas {
49
50// set the default energy scale
51 double GeV = 1.0; //1000;
52
53JetConeFinderTool::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
60JetConeFinderTool::~JetConeFinderTool()
61{}
62
63/////////////////////////////////////////////////////////////////////////////////
64//Execution /
65/////////////////////////////////////////////////////////////////////////////////
66int 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
102void
103JetConeFinderTool::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
235Jet* 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
259FASTJET_END_NAMESPACE
Note: See TracBrowser for help on using the repository browser.