Fork me on GitHub

source: svn/trunk/Utilities/Fastjet/plugins/CDFCones/src/JetCluAlgorithm.cc@ 38

Last change on this file since 38 was 11, checked in by severine ovyn, 16 years ago

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