Fork me on GitHub

source: svn/trunk/Utilities/Fastjet/plugins/CDFCones/src/MidPointAlgorithm.cc@ 14

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

Fastjet added; CDFCones directory has been changed
File size: 12.2 KB
Line 
1#include "../interface/MidPointAlgorithm.hh"
2#include "../interface/ClusterComparisons.hh"
3#include <algorithm>
4#include <iostream>
5#include <cmath>
6
7void MidPointAlgorithm::findStableConesFromSeeds(std::vector<PhysicsTower>& towers, std::vector<Cluster>& stableCones)
8{
9 bool reduceConeSize = true;
10 for(std::vector<PhysicsTower>::iterator towerIter = towers.begin(); towerIter != towers.end(); towerIter++)
11 if(towerIter->fourVector.pt() > _seedThreshold)
12 iterateCone(towerIter->fourVector.y(),towerIter->fourVector.phi(),0,towers,stableCones,reduceConeSize);
13}
14
15void MidPointAlgorithm::findStableConesFromMidPoints(std::vector<PhysicsTower>& towers, std::vector<Cluster>& stableCones)
16{
17 // distanceOK[i-1][j] = Is distance between stableCones i and j (i>j) less than 2*_coneRadius?
18 std::vector< std::vector<bool> > distanceOK;
19 distanceOK.resize(stableCones.size() - 1);
20 for(int nCluster1 = 1; nCluster1 < stableCones.size(); nCluster1++){
21 distanceOK[nCluster1 - 1].resize(nCluster1);
22 double cluster1Rapidity = stableCones[nCluster1].fourVector.y();
23 double cluster1Phi = stableCones[nCluster1].fourVector.phi();
24 for(int nCluster2 = 0; nCluster2 < nCluster1; nCluster2++){
25 double cluster2Rapidity = stableCones[nCluster2].fourVector.y();
26 double cluster2Phi = stableCones[nCluster2].fourVector.phi();
27 double dRapidity = fabs(cluster1Rapidity - cluster2Rapidity);
28 double dPhi = fabs(cluster1Phi - cluster2Phi);
29 if(dPhi > M_PI)
30 dPhi = 2*M_PI - dPhi;
31 double dR = sqrt(dRapidity*dRapidity + dPhi*dPhi);
32 distanceOK[nCluster1 - 1][nCluster2] = dR < 2*_coneRadius;
33 }
34 }
35
36 // Find all pairs (triplets, ...) of stableCones which are less than 2*_coneRadius apart from each other.
37 std::vector< std::vector<int> > pairs(0);
38 std::vector<int> testPair(0);
39 int maxClustersInPair = _maxPairSize;
40 if(!maxClustersInPair)
41 maxClustersInPair = stableCones.size();
42 addClustersToPairs(testPair,pairs,distanceOK,maxClustersInPair);
43
44 // Loop over all combinations. Calculate MidPoint. Make midPointClusters.
45 bool reduceConeSize = false;
46 for(int iPair = 0; iPair < pairs.size(); iPair++){
47 // Calculate rapidity, phi and pT of MidPoint.
48 LorentzVector midPoint(0,0,0,0);
49 for(int iPairMember = 0; iPairMember < pairs[iPair].size(); iPairMember++)
50 midPoint.add(stableCones[pairs[iPair][iPairMember]].fourVector);
51 iterateCone(midPoint.y(),midPoint.phi(),midPoint.pt(),towers,stableCones,reduceConeSize);
52 }
53
54 //sort(stableCones.begin(),stableCones.end(),ClusterPtGreater());
55 local_sort(stableCones); // GPS mod. to allow split-merge with various scales
56}
57
58
59void MidPointAlgorithm::iterateCone(volatile double startRapidity, volatile double startPhi, volatile double startPt,
60 std::vector<PhysicsTower>& towers, std::vector<Cluster>& stableCones, bool reduceConeSize)
61{
62 int nIterations = 0;
63 bool keepJet = true;
64 Cluster trialCone;
65 double iterationConeRadius = _coneRadius;
66 if(reduceConeSize)
67 iterationConeRadius *= sqrt(_coneAreaFraction);
68 while(nIterations++ < _maxIterations + 1 && keepJet){
69 trialCone.clear();
70 // Find particles which should go in the cone.
71 if(nIterations == _maxIterations + 1)
72 iterationConeRadius = _coneRadius;
73 for(std::vector<PhysicsTower>::iterator towerIter = towers.begin(); towerIter != towers.end(); towerIter++){
74 double dRapidity = fabs(towerIter->fourVector.y() - startRapidity);
75 double dPhi = fabs(towerIter->fourVector.phi() - startPhi);
76 if(dPhi > M_PI)
77 dPhi = 2*M_PI - dPhi;
78 double dR = sqrt(dRapidity*dRapidity + dPhi*dPhi);
79 if(dR < iterationConeRadius)
80 trialCone.addTower(*towerIter);
81 }
82 if(!trialCone.size()) // Empty cone?
83 keepJet = false;
84 else{
85 if(nIterations <= _maxIterations){
86 volatile double endRapidity = trialCone.fourVector.y();
87 volatile double endPhi = trialCone.fourVector.phi();
88 volatile double endPt = trialCone.fourVector.pt();
89 // Do we have a stable cone?
90 if(endRapidity == startRapidity && endPhi == startPhi && endPt == startPt){
91 // If cone size is reduced, then do one more iteration.
92 nIterations = _maxIterations;
93 if(!reduceConeSize)
94 nIterations++;
95 }
96 else{
97 // Another iteration.
98 startRapidity = endRapidity;
99 startPhi = endPhi;
100 startPt = endPt;
101 }
102 }
103 }
104 }
105
106 if(keepJet){
107 // We have a stable cone.
108 bool identical = false;
109 for(std::vector<Cluster>::iterator stableConeIter = stableCones.begin(); stableConeIter != stableCones.end(); stableConeIter++)
110 if(trialCone.fourVector.isEqual(stableConeIter->fourVector))
111 identical = true;
112 if(!identical)
113 stableCones.push_back(trialCone);
114 }
115}
116
117void MidPointAlgorithm::addClustersToPairs(std::vector<int>& testPair, std::vector< std::vector<int> >& pairs,
118 std::vector< std::vector<bool> >& distanceOK, int maxClustersInPair)
119{
120 // Recursively adds clusters to pairs, triplets, ... whose mid-points are then calculated.
121
122 // Find StableCone number to start with (either 0 at the beginning or last element of testPair + 1).
123 int nextClusterStart = 0;
124 if(testPair.size())
125 nextClusterStart = testPair.back() + 1;
126 for(int nextCluster = nextClusterStart; nextCluster <= distanceOK.size(); nextCluster++){
127 // Is new SeedCone less than 2*_coneRadius apart from all clusters in testPair?
128 bool addCluster = true;
129 for(int iCluster = 0; iCluster < testPair.size() && addCluster; iCluster++)
130 if(!distanceOK[nextCluster - 1][testPair[iCluster]])
131 addCluster = false;
132 if(addCluster){
133 // Add it to the testPair.
134 testPair.push_back(nextCluster);
135 // If testPair is a pair, add it to pairs.
136 if(testPair.size() > 1)
137 pairs.push_back(testPair);
138 // If not bigger than allowed, find more clusters within 2*_coneRadius.
139 if(testPair.size() < maxClustersInPair)
140 addClustersToPairs(testPair,pairs,distanceOK,maxClustersInPair);
141 // All combinations containing testPair found. Remove last element.
142 testPair.pop_back();
143 }
144 }
145}
146
147void MidPointAlgorithm::splitAndMerge(std::vector<Cluster>& stableCones, std::vector<Cluster>& jets)
148{
149 bool mergingNotFinished = true;
150 while(mergingNotFinished){
151 // Sort the stable cones (highest pt first).
152 //sort(stableCones.begin(),stableCones.end(),ClusterPtGreater());
153 local_sort(stableCones);
154
155 // Start with the highest pt cone.
156 std::vector<Cluster>::iterator stableConeIter1 = stableCones.begin();
157 if(stableConeIter1 == stableCones.end()) // Stable cone list empty?
158 mergingNotFinished = false;
159 else{
160 bool coneNotModified = true;
161 // Determine whether highest pt cone has an overlap with other stable cones.
162 std::vector<Cluster>::iterator stableConeIter2 = stableConeIter1;
163 stableConeIter2++; // 2nd highest pt cone.
164 while(coneNotModified && stableConeIter2 != stableCones.end()){
165 // Calculate overlap of the two cones.
166 Cluster overlap;
167 for(std::vector<PhysicsTower>::iterator towerIter1 = stableConeIter1->towerList.begin();
168 towerIter1 != stableConeIter1->towerList.end();
169 towerIter1++){
170 bool isInCone2 = false;
171 for(std::vector<PhysicsTower>::iterator towerIter2 = stableConeIter2->towerList.begin();
172 towerIter2 != stableConeIter2->towerList.end();
173 towerIter2++)
174 if(towerIter1->isEqual(*towerIter2))
175 isInCone2 = true;
176 if(isInCone2)
177 overlap.addTower(*towerIter1);
178 }
179 if(overlap.size()){ // non-empty overlap
180 coneNotModified = false;
181 // GPS mod to allow various scale choices in split merge --------
182 double overlap_scale, jet2_scale;
183 switch(_smScale) {
184 case SM_pt:
185 overlap_scale = overlap.fourVector.pt();
186 jet2_scale = stableConeIter2->fourVector.pt();
187 break;
188 case SM_Et:
189 overlap_scale = overlap.fourVector.Et();
190 jet2_scale = stableConeIter2->fourVector.Et();
191 break;
192 case SM_mt:
193 overlap_scale = overlap.fourVector.mt();
194 jet2_scale = stableConeIter2->fourVector.mt();
195 break;
196 case SM_pttilde:
197 overlap_scale = overlap.pt_tilde;
198 jet2_scale = stableConeIter2->pt_tilde;
199 break;
200 default:
201 std::cerr << "Unrecognized value for _smScale: "
202 << _smScale << std::endl;
203 exit(-1);
204 }
205 if(overlap_scale >= _overlapThreshold*jet2_scale){
206 // end of GPS modification ---------------------------
207 //if(overlap.fourVector.pt() >= _overlapThreshold*stableConeIter2->fourVector.pt()){
208 // Merge the two cones.
209 for(std::vector<PhysicsTower>::iterator towerIter2 = stableConeIter2->towerList.begin();
210 towerIter2 != stableConeIter2->towerList.end();
211 towerIter2++){
212 bool isInOverlap = false;
213 for(std::vector<PhysicsTower>::iterator overlapTowerIter = overlap.towerList.begin();
214 overlapTowerIter != overlap.towerList.end();
215 overlapTowerIter++)
216 if(towerIter2->isEqual(*overlapTowerIter))
217 isInOverlap = true;
218 if(!isInOverlap)
219 stableConeIter1->addTower(*towerIter2);
220 }
221 // Remove the second cone.
222 stableCones.erase(stableConeIter2);
223 }
224 else{
225 // Separate the two cones.
226 // Which particle goes where?
227 std::vector<PhysicsTower> removeFromCone1,removeFromCone2;
228 for(std::vector<PhysicsTower>::iterator towerIter = overlap.towerList.begin();
229 towerIter != overlap.towerList.end();
230 towerIter++){
231 double towerRapidity = towerIter->fourVector.y();
232 double towerPhi = towerIter->fourVector.phi();
233 // Calculate distance from cone 1.
234 double dRapidity1 = fabs(towerRapidity - stableConeIter1->fourVector.y());
235 double dPhi1 = fabs(towerPhi - stableConeIter1->fourVector.phi());
236 if(dPhi1 > M_PI)
237 dPhi1 = 2*M_PI - dPhi1;
238 double dRJet1 = sqrt(dRapidity1*dRapidity1 + dPhi1*dPhi1);
239 // Calculate distance from cone 2.
240 double dRapidity2 = fabs(towerRapidity - stableConeIter2->fourVector.y());
241 double dPhi2 = fabs(towerPhi - stableConeIter2->fourVector.phi());
242 if(dPhi2 > M_PI)
243 dPhi2 = 2*M_PI - dPhi2;
244 double dRJet2 = sqrt(dRapidity2*dRapidity2 + dPhi2*dPhi2);
245 if(dRJet1 < dRJet2)
246 // Particle is closer to cone 1. To be removed from cone 2.
247 removeFromCone2.push_back(*towerIter);
248 else
249 // Particle is closer to cone 2. To be removed from cone 1.
250 removeFromCone1.push_back(*towerIter);
251 }
252 // Remove particles in the overlap region from the cones to which they have the larger distance.
253 for(std::vector<PhysicsTower>::iterator towerIter = removeFromCone1.begin();
254 towerIter != removeFromCone1.end();
255 towerIter++)
256 stableConeIter1->removeTower(*towerIter);
257 for(std::vector<PhysicsTower>::iterator towerIter = removeFromCone2.begin();
258 towerIter != removeFromCone2.end();
259 towerIter++)
260 stableConeIter2->removeTower(*towerIter);
261 }
262 }
263 stableConeIter2++;
264 }
265 if(coneNotModified){
266 // Cone 1 has no overlap with any of the other cones and can become a jet.
267 jets.push_back(*stableConeIter1);
268 stableCones.erase(stableConeIter1);
269 }
270 }
271 }
272
273 //sort(jets.begin(),jets.end(),ClusterPtGreater());
274 local_sort(jets); // GPS mod. to allow split-merge with various scales
275}
276
277
278
279void MidPointAlgorithm::local_sort(std::vector<Cluster>& clusters) {
280 switch(_smScale) {
281 case SM_pt:
282 sort(clusters.begin(),clusters.end(),ClusterPtGreater());
283 break;
284 case SM_Et:
285 sort(clusters.begin(),clusters.end(),ClusterFourVectorEtGreater());
286 break;
287 case SM_mt:
288 sort(clusters.begin(),clusters.end(),ClusterMtGreater());
289 break;
290 case SM_pttilde:
291 sort(clusters.begin(),clusters.end(),ClusterPtTildeGreater());
292 break;
293 default:
294 std::cerr << "Unrecognized value for _smScale: " << _smScale << std::endl;
295 exit(-1);
296 }
297}
298
299
300
301void MidPointAlgorithm::run(std::vector<PhysicsTower>& towers, std::vector<Cluster>& jets)
302{
303 std::vector<Cluster> stableCones;
304 findStableConesFromSeeds(towers,stableCones);
305 // GPS addition to prevent crashes if no stable cones
306 // are found (e.g. all particles below seed threshold)
307 if (stableCones.size() > 0) {
308 findStableConesFromMidPoints(towers,stableCones);
309 splitAndMerge(stableCones,jets);
310 }
311}
312
313
Note: See TracBrowser for help on using the repository browser.