Fork me on GitHub

source: svn/trunk/external/fastjet/plugins/CDFCones/JetCluAlgorithm.cc@ 1250

Last change on this file since 1250 was 859, checked in by Pavel Demin, 12 years ago

update fastjet to version 3.0.3
  • Property svn:eol-style set to native
  • Property svn:keywords set to Id Revision Date
File size: 13.3 KB
Line 
1//----------------------------------------------------------------------
2// This file distributed with FastJet has been obtained from
3// http://www.pa.msu.edu/~huston/Les_Houches_2005/JetClu+Midpoint-StandAlone.tgz
4//
5// Permission to distribute it with FastJet has been granted by Joey
6// Huston (see the COPYING file in the main FastJet directory for
7// details).
8// Changes from the original file are listed below.
9//----------------------------------------------------------------------
10
11// History of changes compared to the original JetCluAlgorithm.cc file
12//
13// 2011-11-14 Gregory Soyez <soyez@fastjet.fr>
14//
15// * added a few parentheses suggested by the -Wparentheses gcc option
16//
17// 2009-01-17 Gregory Soyez <soyez@fastjet.fr>
18//
19// * put the code in the fastjet::cdf namespace
20//
21// 2006-09-24 Gavin Salam <salam@lpthe.jussieu.fr>
22//
23// * added JetClu+MidPoint to FastJet
24
25#include "JetCluAlgorithm.hh"
26#include "ClusterComparisons.hh"
27#include "Centroid.hh"
28#include <algorithm>
29#include <cmath>
30
31#include <fastjet/internal/base.hh>
32
33FASTJET_BEGIN_NAMESPACE
34
35namespace cdf{
36
37void JetCluAlgorithm::makeSeedTowers(std::vector<PhysicsTower>& towers, std::vector<Cluster>& seedTowers)
38{
39 for(int iEta = 4; iEta < 48; iEta++){
40 bool seg24 = true;
41 if ((iEta >= 8 && iEta < 14) || (iEta >= 38 && iEta < 44))
42 seg24 = false;
43 for(int iPhi = 0; iPhi < 24; iPhi++){
44 Cluster seed;
45 for(std::vector<PhysicsTower>::iterator towerIter = towers.begin(); towerIter != towers.end(); towerIter++)
46 if(towerIter->iEta() == iEta &&
47 ((seg24 && towerIter->iPhi() == iPhi) || (!seg24 && (towerIter->iPhi() == 2*iPhi || towerIter->iPhi() == 2*iPhi + 1))))
48 seed.addTower(*towerIter);
49 if(seed.centroid.Et > _seedThreshold)
50 seedTowers.push_back(seed);
51 }
52 }
53 sort(seedTowers.begin(),seedTowers.end(),ClusterCentroidEtGreater());
54}
55
56void JetCluAlgorithm::buildPreClusters(std::vector<Cluster>& seedTowers, std::vector<PhysicsTower>& towers,
57 std::vector<Cluster>& preClusters)
58{
59 std::vector<Centroid> leadingSeedTowers;
60 for(std::vector<Cluster>::iterator seedTowerIter = seedTowers.begin(); seedTowerIter != seedTowers.end(); seedTowerIter++){
61 bool seedTowerAddedToPreCluster = false;
62 std::vector<Cluster>::iterator preClusterIter = preClusters.begin();
63 std::vector<Centroid>::iterator leadingSeedTowerIter = leadingSeedTowers.begin();
64 while(preClusterIter != preClusters.end() && !seedTowerAddedToPreCluster){
65 double dEta = fabs(seedTowerIter->centroid.eta - leadingSeedTowerIter->eta);
66 double dPhi = fabs(seedTowerIter->centroid.phi - leadingSeedTowerIter->phi);
67 if(dPhi > M_PI)
68 dPhi = 2*M_PI - dPhi;
69 if(dEta <= _coneRadius && dPhi <= _coneRadius){
70 int iEtaSeedTower = seedTowerIter->towerList.begin()->iEta();
71 int iPhiSeedTower = seedTowerIter->towerList.begin()->iPhi();
72 if ((iEtaSeedTower >= 8 && iEtaSeedTower < 14) || (iEtaSeedTower >= 38 && iEtaSeedTower < 44))
73 iPhiSeedTower = iPhiSeedTower/2;
74 for(std::vector<PhysicsTower>::iterator preClusterTowerIter = preClusterIter->towerList.begin();
75 preClusterTowerIter != preClusterIter->towerList.end() && !seedTowerAddedToPreCluster;
76 preClusterTowerIter++){
77 int iEtaPreClusterTower = preClusterTowerIter->iEta();
78 int iPhiPreClusterTower = preClusterTowerIter->iPhi();
79 if ((iEtaPreClusterTower >= 8 && iEtaPreClusterTower < 14) || (iEtaPreClusterTower >= 38 && iEtaPreClusterTower < 44))
80 iPhiPreClusterTower = iPhiPreClusterTower/2;
81 int dIEta = abs(iEtaSeedTower - iEtaPreClusterTower);
82 int dIPhi = abs(iPhiSeedTower - iPhiPreClusterTower);
83 if(dIPhi > 12)
84 dIPhi = 24 - dIPhi;
85 int adj = dIPhi*dIPhi + dIEta*dIEta;
86 if(adj <= _adjacencyCut){
87 for(std::vector<PhysicsTower>::iterator seedTowerTowerIter = seedTowerIter->towerList.begin();
88 seedTowerTowerIter != seedTowerIter->towerList.end();
89 seedTowerTowerIter++)
90 preClusterIter->addTower(*seedTowerTowerIter);
91 seedTowerAddedToPreCluster = true;
92 }
93 }
94 }
95 preClusterIter++;
96 leadingSeedTowerIter++;
97 }
98 if(!seedTowerAddedToPreCluster){
99 Cluster newPreCluster;
100 for(std::vector<PhysicsTower>::iterator seedTowerTowerIter = seedTowerIter->towerList.begin();
101 seedTowerTowerIter != seedTowerIter->towerList.end();
102 seedTowerTowerIter++)
103 newPreCluster.addTower(*seedTowerTowerIter);
104 preClusters.push_back(newPreCluster);
105 leadingSeedTowers.push_back(Centroid(newPreCluster.centroid.Et,newPreCluster.centroid.eta,newPreCluster.centroid.phi));
106 }
107 }
108}
109
110void JetCluAlgorithm::findStableCones(std::vector<Cluster>& preClusters, std::vector<PhysicsTower>& towers,
111 std::vector<Cluster>& stableCones)
112{
113 for(std::vector<Cluster>::iterator preClusterIter = preClusters.begin(); preClusterIter != preClusters.end(); preClusterIter++){
114 double startEt = preClusterIter->centroid.Et;
115 double startEta = preClusterIter->centroid.eta;
116 double startPhi = preClusterIter->centroid.phi;
117 int nIterations = 0;
118 Cluster trialCone;
119 while(nIterations++ < _maxIterations){
120 trialCone.clear();
121 for(std::vector<PhysicsTower>::iterator towerIter = towers.begin(); towerIter != towers.end(); towerIter++){
122 double dEta = fabs(towerIter->eta() - startEta);
123 double dPhi = fabs(towerIter->phi() - startPhi);
124 if(dPhi > M_PI)
125 dPhi = 2*M_PI - dPhi;
126 double dR = sqrt(dEta*dEta + dPhi*dPhi);
127 if(dR < _coneRadius)
128 trialCone.addTower(*towerIter);
129 }
130 if(_iratch != 0)
131 for(std::vector<PhysicsTower>::iterator preClusterTowerIter = preClusterIter->towerList.begin();
132 preClusterTowerIter != preClusterIter->towerList.end();
133 preClusterTowerIter++){
134 bool foundInTrialCone = false;
135 for(std::vector<PhysicsTower>::iterator trialConeTowerIter = trialCone.towerList.begin();
136 trialConeTowerIter != trialCone.towerList.end() && !foundInTrialCone;
137 trialConeTowerIter++)
138 if(trialConeTowerIter->isEqual(*preClusterTowerIter))
139 foundInTrialCone = true;
140 if(!foundInTrialCone)
141 trialCone.addTower(*preClusterTowerIter);
142 }
143 if(nIterations <= _maxIterations){
144 double endEt = trialCone.centroid.Et;
145 double endEta = trialCone.centroid.eta;
146 double endPhi = trialCone.centroid.phi;
147 if(endEt == startEt && endEta == startEta && endPhi == startPhi)
148 nIterations = _maxIterations;
149 else{
150 startEt = endEt;
151 startEta = endEta;
152 startPhi = endPhi;
153 }
154 }
155 }
156// bool foundIdentical = false;
157// for(std::vector<Cluster>::iterator stableConeIter = stableCones.begin();
158// stableConeIter != stableCones.end() && !foundIdentical;
159// stableConeIter++)
160// if(trialCone.centroid.isEqual(stableConeIter->centroid))
161// foundIdentical = true;
162// if(!foundIdentical)
163 stableCones.push_back(trialCone);
164 }
165 sort(stableCones.begin(),stableCones.end(),ClusterCentroidEtGreater());
166}
167
168void JetCluAlgorithm::splitAndMerge(std::vector<Cluster>& stableCones, std::vector<Cluster>& jets)
169{
170 std::vector<bool> isActive;
171 for(std::vector<Cluster>::iterator stableConeIter = stableCones.begin(); stableConeIter != stableCones.end(); stableConeIter++)
172 isActive.push_back(bool(true));
173 std::vector<bool>::iterator isActiveIter1 = isActive.begin();
174 for(std::vector<Cluster>::iterator stableConeIter1 = stableCones.begin();
175 stableConeIter1 != stableCones.end();
176 stableConeIter1++, isActiveIter1++){
177 std::vector<Cluster>::iterator stableConeIter2 = stableCones.begin();
178 std::vector<bool>::iterator isActiveIter2 = isActive.begin();
179 while(stableConeIter2 != stableConeIter1 && *isActiveIter1){
180 if(*isActiveIter2){
181 Cluster overlap;
182 for(std::vector<PhysicsTower>::iterator towerIter1 = stableConeIter1->towerList.begin();
183 towerIter1 != stableConeIter1->towerList.end();
184 towerIter1++)
185 for(std::vector<PhysicsTower>::iterator towerIter2 = stableConeIter2->towerList.begin();
186 towerIter2 != stableConeIter2->towerList.end();
187 towerIter2++)
188 if(towerIter1->isEqual(*towerIter2)){
189 overlap.addTower(*towerIter1);
190 break;
191 }
192 if(overlap.size()){
193 if(overlap.size() == stableConeIter2->size())
194 *isActiveIter2 = false;
195 else if(overlap.size() == stableConeIter1->size())
196 *isActiveIter1 = false;
197 else if(overlap.centroid.Et > _overlapThreshold*stableConeIter1->centroid.Et ||
198 overlap.centroid.Et > _overlapThreshold*stableConeIter2->centroid.Et){
199 for(std::vector<PhysicsTower>::iterator stableConeTowerIter2 = stableConeIter2->towerList.begin();
200 stableConeTowerIter2 != stableConeIter2->towerList.end();
201 stableConeTowerIter2++){
202 bool isInOverlap = false;
203 for(std::vector<PhysicsTower>::iterator overlapTowerIter = overlap.towerList.begin();
204 overlapTowerIter != overlap.towerList.end() && !isInOverlap;
205 overlapTowerIter++)
206 if(stableConeTowerIter2->isEqual(*overlapTowerIter))
207 isInOverlap = true;
208 if(!isInOverlap)
209 stableConeIter1->addTower(*stableConeTowerIter2);
210 }
211 *isActiveIter2 = false;
212 }
213 else{
214 Cluster removeFromStableCone1,removeFromStableCone2,oldRemoveFromStableCone1,oldRemoveFromStableCone2;
215 double etaStableCone1 = stableConeIter1->centroid.eta;
216 double phiStableCone1 = stableConeIter1->centroid.phi;
217 double etaStableCone2 = stableConeIter2->centroid.eta;
218 double phiStableCone2 = stableConeIter2->centroid.phi;
219 double dRstableCone1,dRstableCone2;
220 int iterCount = 0;
221 while(iterCount++ <= _maxIterations){
222 oldRemoveFromStableCone1.clear();
223 oldRemoveFromStableCone2.clear();
224 if(iterCount > 1){
225 if(removeFromStableCone1.size()){
226 Centroid stableConeCentroid1(stableConeIter1->centroid);
227 Centroid removeCentroid1(removeFromStableCone1.centroid);
228 stableConeCentroid1.subtract(removeCentroid1);
229 etaStableCone1 = stableConeCentroid1.eta;
230 phiStableCone1 = stableConeCentroid1.phi;
231 }
232 else{
233 etaStableCone1 = stableConeIter1->centroid.eta;
234 phiStableCone1 = stableConeIter1->centroid.phi;
235 }
236 if(removeFromStableCone2.size()){
237 Centroid stableConeCentroid2(stableConeIter2->centroid);
238 Centroid removeCentroid2(removeFromStableCone2.centroid);
239 stableConeCentroid2.subtract(removeCentroid2);
240 etaStableCone2 = stableConeCentroid2.eta;
241 phiStableCone2 = stableConeCentroid2.phi;
242 }
243 else{
244 etaStableCone2 = stableConeIter2->centroid.eta;
245 phiStableCone2 = stableConeIter2->centroid.phi;
246 }
247 for(std::vector<PhysicsTower>::iterator removeTowerIter1 = removeFromStableCone1.towerList.begin();
248 removeTowerIter1 != removeFromStableCone1.towerList.end();
249 removeTowerIter1++)
250 oldRemoveFromStableCone1.addTower(*removeTowerIter1);
251 for(std::vector<PhysicsTower>::iterator removeTowerIter2 = removeFromStableCone2.towerList.begin();
252 removeTowerIter2 != removeFromStableCone2.towerList.end();
253 removeTowerIter2++)
254 oldRemoveFromStableCone2.addTower(*removeTowerIter2);
255 }
256 removeFromStableCone1.clear();
257 removeFromStableCone2.clear();
258 for(std::vector<PhysicsTower>::iterator overlapTowerIter = overlap.towerList.begin();
259 overlapTowerIter != overlap.towerList.end();
260 overlapTowerIter++){
261 double dEta1 = fabs(overlapTowerIter->eta() - etaStableCone1);
262 double dPhi1 = fabs(overlapTowerIter->phi() - phiStableCone1);
263 if(dPhi1 > M_PI)
264 dPhi1 = 2*M_PI - dPhi1;
265 dRstableCone1 = dEta1*dEta1 + dPhi1*dPhi1;
266 double dEta2 = fabs(overlapTowerIter->eta() - etaStableCone2);
267 double dPhi2 = fabs(overlapTowerIter->phi() - phiStableCone2);
268 if(dPhi2 > M_PI)
269 dPhi2 = 2*M_PI - dPhi2;
270 dRstableCone2 = dEta2*dEta2 + dPhi2*dPhi2;
271 if(dRstableCone1 < dRstableCone2)
272 removeFromStableCone2.addTower(*overlapTowerIter);
273 else
274 removeFromStableCone1.addTower(*overlapTowerIter);
275 }
276 if(iterCount > 1 &&
277 removeFromStableCone1.size() == oldRemoveFromStableCone1.size() &&
278 removeFromStableCone2.size() == oldRemoveFromStableCone2.size() &&
279 (!removeFromStableCone1.size() || !removeFromStableCone2.size() ||
280 (removeFromStableCone1.centroid.isEqual(oldRemoveFromStableCone1.centroid) &&
281 removeFromStableCone2.centroid.isEqual(oldRemoveFromStableCone2.centroid))))
282 iterCount = _maxIterations + 1;
283 }
284 for(std::vector<PhysicsTower>::iterator removeTowerIter1 = removeFromStableCone1.towerList.begin();
285 removeTowerIter1 != removeFromStableCone1.towerList.end();
286 removeTowerIter1++)
287 stableConeIter1->removeTower(*removeTowerIter1);
288 for(std::vector<PhysicsTower>::iterator removeTowerIter2 = removeFromStableCone2.towerList.begin();
289 removeTowerIter2 != removeFromStableCone2.towerList.end();
290 removeTowerIter2++)
291 stableConeIter2->removeTower(*removeTowerIter2);
292 }
293 overlap.clear();
294 }
295 }
296 stableConeIter2++;
297 isActiveIter2++;
298 }
299 }
300 jets.clear();
301 std::vector<bool>::iterator isActiveIter = isActive.begin();
302 for(std::vector<Cluster>::iterator stableConeIter = stableCones.begin();
303 stableConeIter != stableCones.end();
304 stableConeIter++, isActiveIter++)
305 if(*isActiveIter)
306 jets.push_back(*stableConeIter);
307 sort(jets.begin(),jets.end(),ClusterFourVectorEtGreater());
308}
309
310void JetCluAlgorithm::run(std::vector<PhysicsTower>& towers, std::vector<Cluster>& jets)
311{
312 std::vector<Cluster> seedTowers;
313 makeSeedTowers(towers,seedTowers);
314 std::vector<Cluster> preClusters;
315 buildPreClusters(seedTowers,towers,preClusters);
316 std::vector<Cluster> stableCones;
317 findStableCones(preClusters,towers,stableCones);
318 splitAndMerge(stableCones,jets);
319}
320
321} // namespace cdf
322
323FASTJET_END_NAMESPACE
Note: See TracBrowser for help on using the repository browser.