Fork me on GitHub

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

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

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