Fork me on GitHub

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

Last change on this file since 1377 was 1332, checked in by Pavel Demin, 11 years ago

upgrade fastjet to 3.0.6
  • 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#include <cstdlib>
31
32#include <fastjet/internal/base.hh>
33
34FASTJET_BEGIN_NAMESPACE
35
36namespace cdf{
37
38void JetCluAlgorithm::makeSeedTowers(std::vector<PhysicsTower>& towers, std::vector<Cluster>& seedTowers)
39{
40 for(int iEta = 4; iEta < 48; iEta++){
41 bool seg24 = true;
42 if ((iEta >= 8 && iEta < 14) || (iEta >= 38 && iEta < 44))
43 seg24 = false;
44 for(int iPhi = 0; iPhi < 24; iPhi++){
45 Cluster seed;
46 for(std::vector<PhysicsTower>::iterator towerIter = towers.begin(); towerIter != towers.end(); towerIter++)
47 if(towerIter->iEta() == iEta &&
48 ((seg24 && towerIter->iPhi() == iPhi) || (!seg24 && (towerIter->iPhi() == 2*iPhi || towerIter->iPhi() == 2*iPhi + 1))))
49 seed.addTower(*towerIter);
50 if(seed.centroid.Et > _seedThreshold)
51 seedTowers.push_back(seed);
52 }
53 }
54 sort(seedTowers.begin(),seedTowers.end(),ClusterCentroidEtGreater());
55}
56
57void JetCluAlgorithm::buildPreClusters(std::vector<Cluster>& seedTowers, std::vector<PhysicsTower>& towers,
58 std::vector<Cluster>& preClusters)
59{
60 std::vector<Centroid> leadingSeedTowers;
61 for(std::vector<Cluster>::iterator seedTowerIter = seedTowers.begin(); seedTowerIter != seedTowers.end(); seedTowerIter++){
62 bool seedTowerAddedToPreCluster = false;
63 std::vector<Cluster>::iterator preClusterIter = preClusters.begin();
64 std::vector<Centroid>::iterator leadingSeedTowerIter = leadingSeedTowers.begin();
65 while(preClusterIter != preClusters.end() && !seedTowerAddedToPreCluster){
66 double dEta = fabs(seedTowerIter->centroid.eta - leadingSeedTowerIter->eta);
67 double dPhi = fabs(seedTowerIter->centroid.phi - leadingSeedTowerIter->phi);
68 if(dPhi > M_PI)
69 dPhi = 2*M_PI - dPhi;
70 if(dEta <= _coneRadius && dPhi <= _coneRadius){
71 int iEtaSeedTower = seedTowerIter->towerList.begin()->iEta();
72 int iPhiSeedTower = seedTowerIter->towerList.begin()->iPhi();
73 if ((iEtaSeedTower >= 8 && iEtaSeedTower < 14) || (iEtaSeedTower >= 38 && iEtaSeedTower < 44))
74 iPhiSeedTower = iPhiSeedTower/2;
75 for(std::vector<PhysicsTower>::iterator preClusterTowerIter = preClusterIter->towerList.begin();
76 preClusterTowerIter != preClusterIter->towerList.end() && !seedTowerAddedToPreCluster;
77 preClusterTowerIter++){
78 int iEtaPreClusterTower = preClusterTowerIter->iEta();
79 int iPhiPreClusterTower = preClusterTowerIter->iPhi();
80 if ((iEtaPreClusterTower >= 8 && iEtaPreClusterTower < 14) || (iEtaPreClusterTower >= 38 && iEtaPreClusterTower < 44))
81 iPhiPreClusterTower = iPhiPreClusterTower/2;
82 int dIEta = std::abs(iEtaSeedTower - iEtaPreClusterTower);
83 int dIPhi = std::abs(iPhiSeedTower - iPhiPreClusterTower);
84 if(dIPhi > 12)
85 dIPhi = 24 - dIPhi;
86 int adj = dIPhi*dIPhi + dIEta*dIEta;
87 if(adj <= _adjacencyCut){
88 for(std::vector<PhysicsTower>::iterator seedTowerTowerIter = seedTowerIter->towerList.begin();
89 seedTowerTowerIter != seedTowerIter->towerList.end();
90 seedTowerTowerIter++)
91 preClusterIter->addTower(*seedTowerTowerIter);
92 seedTowerAddedToPreCluster = true;
93 }
94 }
95 }
96 preClusterIter++;
97 leadingSeedTowerIter++;
98 }
99 if(!seedTowerAddedToPreCluster){
100 Cluster newPreCluster;
101 for(std::vector<PhysicsTower>::iterator seedTowerTowerIter = seedTowerIter->towerList.begin();
102 seedTowerTowerIter != seedTowerIter->towerList.end();
103 seedTowerTowerIter++)
104 newPreCluster.addTower(*seedTowerTowerIter);
105 preClusters.push_back(newPreCluster);
106 leadingSeedTowers.push_back(Centroid(newPreCluster.centroid.Et,newPreCluster.centroid.eta,newPreCluster.centroid.phi));
107 }
108 }
109}
110
111void JetCluAlgorithm::findStableCones(std::vector<Cluster>& preClusters, std::vector<PhysicsTower>& towers,
112 std::vector<Cluster>& stableCones)
113{
114 for(std::vector<Cluster>::iterator preClusterIter = preClusters.begin(); preClusterIter != preClusters.end(); preClusterIter++){
115 double startEt = preClusterIter->centroid.Et;
116 double startEta = preClusterIter->centroid.eta;
117 double startPhi = preClusterIter->centroid.phi;
118 int nIterations = 0;
119 Cluster trialCone;
120 while(nIterations++ < _maxIterations){
121 trialCone.clear();
122 for(std::vector<PhysicsTower>::iterator towerIter = towers.begin(); towerIter != towers.end(); towerIter++){
123 double dEta = fabs(towerIter->eta() - startEta);
124 double dPhi = fabs(towerIter->phi() - startPhi);
125 if(dPhi > M_PI)
126 dPhi = 2*M_PI - dPhi;
127 double dR = sqrt(dEta*dEta + dPhi*dPhi);
128 if(dR < _coneRadius)
129 trialCone.addTower(*towerIter);
130 }
131 if(_iratch != 0)
132 for(std::vector<PhysicsTower>::iterator preClusterTowerIter = preClusterIter->towerList.begin();
133 preClusterTowerIter != preClusterIter->towerList.end();
134 preClusterTowerIter++){
135 bool foundInTrialCone = false;
136 for(std::vector<PhysicsTower>::iterator trialConeTowerIter = trialCone.towerList.begin();
137 trialConeTowerIter != trialCone.towerList.end() && !foundInTrialCone;
138 trialConeTowerIter++)
139 if(trialConeTowerIter->isEqual(*preClusterTowerIter))
140 foundInTrialCone = true;
141 if(!foundInTrialCone)
142 trialCone.addTower(*preClusterTowerIter);
143 }
144 if(nIterations <= _maxIterations){
145 double endEt = trialCone.centroid.Et;
146 double endEta = trialCone.centroid.eta;
147 double endPhi = trialCone.centroid.phi;
148 if(endEt == startEt && endEta == startEta && endPhi == startPhi)
149 nIterations = _maxIterations;
150 else{
151 startEt = endEt;
152 startEta = endEta;
153 startPhi = endPhi;
154 }
155 }
156 }
157// bool foundIdentical = false;
158// for(std::vector<Cluster>::iterator stableConeIter = stableCones.begin();
159// stableConeIter != stableCones.end() && !foundIdentical;
160// stableConeIter++)
161// if(trialCone.centroid.isEqual(stableConeIter->centroid))
162// foundIdentical = true;
163// if(!foundIdentical)
164 stableCones.push_back(trialCone);
165 }
166 sort(stableCones.begin(),stableCones.end(),ClusterCentroidEtGreater());
167}
168
169void JetCluAlgorithm::splitAndMerge(std::vector<Cluster>& stableCones, std::vector<Cluster>& jets)
170{
171 std::vector<bool> isActive;
172 for(std::vector<Cluster>::iterator stableConeIter = stableCones.begin(); stableConeIter != stableCones.end(); stableConeIter++)
173 isActive.push_back(bool(true));
174 std::vector<bool>::iterator isActiveIter1 = isActive.begin();
175 for(std::vector<Cluster>::iterator stableConeIter1 = stableCones.begin();
176 stableConeIter1 != stableCones.end();
177 stableConeIter1++, isActiveIter1++){
178 std::vector<Cluster>::iterator stableConeIter2 = stableCones.begin();
179 std::vector<bool>::iterator isActiveIter2 = isActive.begin();
180 while(stableConeIter2 != stableConeIter1 && *isActiveIter1){
181 if(*isActiveIter2){
182 Cluster overlap;
183 for(std::vector<PhysicsTower>::iterator towerIter1 = stableConeIter1->towerList.begin();
184 towerIter1 != stableConeIter1->towerList.end();
185 towerIter1++)
186 for(std::vector<PhysicsTower>::iterator towerIter2 = stableConeIter2->towerList.begin();
187 towerIter2 != stableConeIter2->towerList.end();
188 towerIter2++)
189 if(towerIter1->isEqual(*towerIter2)){
190 overlap.addTower(*towerIter1);
191 break;
192 }
193 if(overlap.size()){
194 if(overlap.size() == stableConeIter2->size())
195 *isActiveIter2 = false;
196 else if(overlap.size() == stableConeIter1->size())
197 *isActiveIter1 = false;
198 else if(overlap.centroid.Et > _overlapThreshold*stableConeIter1->centroid.Et ||
199 overlap.centroid.Et > _overlapThreshold*stableConeIter2->centroid.Et){
200 for(std::vector<PhysicsTower>::iterator stableConeTowerIter2 = stableConeIter2->towerList.begin();
201 stableConeTowerIter2 != stableConeIter2->towerList.end();
202 stableConeTowerIter2++){
203 bool isInOverlap = false;
204 for(std::vector<PhysicsTower>::iterator overlapTowerIter = overlap.towerList.begin();
205 overlapTowerIter != overlap.towerList.end() && !isInOverlap;
206 overlapTowerIter++)
207 if(stableConeTowerIter2->isEqual(*overlapTowerIter))
208 isInOverlap = true;
209 if(!isInOverlap)
210 stableConeIter1->addTower(*stableConeTowerIter2);
211 }
212 *isActiveIter2 = false;
213 }
214 else{
215 Cluster removeFromStableCone1,removeFromStableCone2,oldRemoveFromStableCone1,oldRemoveFromStableCone2;
216 double etaStableCone1 = stableConeIter1->centroid.eta;
217 double phiStableCone1 = stableConeIter1->centroid.phi;
218 double etaStableCone2 = stableConeIter2->centroid.eta;
219 double phiStableCone2 = stableConeIter2->centroid.phi;
220 double dRstableCone1,dRstableCone2;
221 int iterCount = 0;
222 while(iterCount++ <= _maxIterations){
223 oldRemoveFromStableCone1.clear();
224 oldRemoveFromStableCone2.clear();
225 if(iterCount > 1){
226 if(removeFromStableCone1.size()){
227 Centroid stableConeCentroid1(stableConeIter1->centroid);
228 Centroid removeCentroid1(removeFromStableCone1.centroid);
229 stableConeCentroid1.subtract(removeCentroid1);
230 etaStableCone1 = stableConeCentroid1.eta;
231 phiStableCone1 = stableConeCentroid1.phi;
232 }
233 else{
234 etaStableCone1 = stableConeIter1->centroid.eta;
235 phiStableCone1 = stableConeIter1->centroid.phi;
236 }
237 if(removeFromStableCone2.size()){
238 Centroid stableConeCentroid2(stableConeIter2->centroid);
239 Centroid removeCentroid2(removeFromStableCone2.centroid);
240 stableConeCentroid2.subtract(removeCentroid2);
241 etaStableCone2 = stableConeCentroid2.eta;
242 phiStableCone2 = stableConeCentroid2.phi;
243 }
244 else{
245 etaStableCone2 = stableConeIter2->centroid.eta;
246 phiStableCone2 = stableConeIter2->centroid.phi;
247 }
248 for(std::vector<PhysicsTower>::iterator removeTowerIter1 = removeFromStableCone1.towerList.begin();
249 removeTowerIter1 != removeFromStableCone1.towerList.end();
250 removeTowerIter1++)
251 oldRemoveFromStableCone1.addTower(*removeTowerIter1);
252 for(std::vector<PhysicsTower>::iterator removeTowerIter2 = removeFromStableCone2.towerList.begin();
253 removeTowerIter2 != removeFromStableCone2.towerList.end();
254 removeTowerIter2++)
255 oldRemoveFromStableCone2.addTower(*removeTowerIter2);
256 }
257 removeFromStableCone1.clear();
258 removeFromStableCone2.clear();
259 for(std::vector<PhysicsTower>::iterator overlapTowerIter = overlap.towerList.begin();
260 overlapTowerIter != overlap.towerList.end();
261 overlapTowerIter++){
262 double dEta1 = fabs(overlapTowerIter->eta() - etaStableCone1);
263 double dPhi1 = fabs(overlapTowerIter->phi() - phiStableCone1);
264 if(dPhi1 > M_PI)
265 dPhi1 = 2*M_PI - dPhi1;
266 dRstableCone1 = dEta1*dEta1 + dPhi1*dPhi1;
267 double dEta2 = fabs(overlapTowerIter->eta() - etaStableCone2);
268 double dPhi2 = fabs(overlapTowerIter->phi() - phiStableCone2);
269 if(dPhi2 > M_PI)
270 dPhi2 = 2*M_PI - dPhi2;
271 dRstableCone2 = dEta2*dEta2 + dPhi2*dPhi2;
272 if(dRstableCone1 < dRstableCone2)
273 removeFromStableCone2.addTower(*overlapTowerIter);
274 else
275 removeFromStableCone1.addTower(*overlapTowerIter);
276 }
277 if(iterCount > 1 &&
278 removeFromStableCone1.size() == oldRemoveFromStableCone1.size() &&
279 removeFromStableCone2.size() == oldRemoveFromStableCone2.size() &&
280 (!removeFromStableCone1.size() || !removeFromStableCone2.size() ||
281 (removeFromStableCone1.centroid.isEqual(oldRemoveFromStableCone1.centroid) &&
282 removeFromStableCone2.centroid.isEqual(oldRemoveFromStableCone2.centroid))))
283 iterCount = _maxIterations + 1;
284 }
285 for(std::vector<PhysicsTower>::iterator removeTowerIter1 = removeFromStableCone1.towerList.begin();
286 removeTowerIter1 != removeFromStableCone1.towerList.end();
287 removeTowerIter1++)
288 stableConeIter1->removeTower(*removeTowerIter1);
289 for(std::vector<PhysicsTower>::iterator removeTowerIter2 = removeFromStableCone2.towerList.begin();
290 removeTowerIter2 != removeFromStableCone2.towerList.end();
291 removeTowerIter2++)
292 stableConeIter2->removeTower(*removeTowerIter2);
293 }
294 overlap.clear();
295 }
296 }
297 stableConeIter2++;
298 isActiveIter2++;
299 }
300 }
301 jets.clear();
302 std::vector<bool>::iterator isActiveIter = isActive.begin();
303 for(std::vector<Cluster>::iterator stableConeIter = stableCones.begin();
304 stableConeIter != stableCones.end();
305 stableConeIter++, isActiveIter++)
306 if(*isActiveIter)
307 jets.push_back(*stableConeIter);
308 sort(jets.begin(),jets.end(),ClusterFourVectorEtGreater());
309}
310
311void JetCluAlgorithm::run(std::vector<PhysicsTower>& towers, std::vector<Cluster>& jets)
312{
313 std::vector<Cluster> seedTowers;
314 makeSeedTowers(towers,seedTowers);
315 std::vector<Cluster> preClusters;
316 buildPreClusters(seedTowers,towers,preClusters);
317 std::vector<Cluster> stableCones;
318 findStableCones(preClusters,towers,stableCones);
319 splitAndMerge(stableCones,jets);
320}
321
322} // namespace cdf
323
324FASTJET_END_NAMESPACE
Note: See TracBrowser for help on using the repository browser.