Fork me on GitHub

source: git/external/fastjet/plugins/CDFCones/JetCluAlgorithm.cc@ 4e8e72b

Last change on this file since 4e8e72b was 35cdc46, checked in by Pavel Demin <demin@…>, 10 years ago

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