1 | #include "Utilities/CDFCones/interface/MidPointAlgorithm.h"
|
---|
2 | #include "Utilities/CDFCones/interface/ClusterComparisons.h"
|
---|
3 | #include <algorithm>
|
---|
4 | #include <cmath>
|
---|
5 |
|
---|
6 | void MidPointAlgorithm::findStableConesFromSeeds(std::vector<PhysicsTower>& towers, std::vector<Cluster>& stableCones)
|
---|
7 | {
|
---|
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 |
|
---|
15 | void MidPointAlgorithm::findStableConesFromMidPoints(std::vector<PhysicsTower>& towers, std::vector<Cluster>& stableCones)
|
---|
16 | {
|
---|
17 |
|
---|
18 | // distanceOK[i-1][j] = Is distance between stableCones i and j (i>j) less than 2*_coneRadius?
|
---|
19 | std::vector< std::vector<bool> > distanceOK;
|
---|
20 | if(stableCones.size()>0){
|
---|
21 | distanceOK.resize(stableCones.size() - 1);
|
---|
22 | for(unsigned int nCluster1 = 1; nCluster1 < stableCones.size(); nCluster1++){
|
---|
23 | distanceOK[nCluster1 - 1].resize(nCluster1);
|
---|
24 | double cluster1Rapidity = stableCones[nCluster1].fourVector.y();
|
---|
25 | double cluster1Phi = stableCones[nCluster1].fourVector.phi();
|
---|
26 | for(unsigned int nCluster2 = 0; nCluster2 < nCluster1; nCluster2++){
|
---|
27 | double cluster2Rapidity = stableCones[nCluster2].fourVector.y();
|
---|
28 | double cluster2Phi = stableCones[nCluster2].fourVector.phi();
|
---|
29 | double dRapidity = fabs(cluster1Rapidity - cluster2Rapidity);
|
---|
30 | double dPhi = fabs(cluster1Phi - cluster2Phi);
|
---|
31 | if(dPhi > M_PI)
|
---|
32 | dPhi = 2*M_PI - dPhi;
|
---|
33 | double dR = sqrt(dRapidity*dRapidity + dPhi*dPhi);
|
---|
34 | distanceOK[nCluster1 - 1][nCluster2] = dR < 2*_coneRadius;
|
---|
35 | }
|
---|
36 | }
|
---|
37 | }
|
---|
38 | // Find all pairs (triplets, ...) of stableCones which are less than 2*_coneRadius apart from each other.
|
---|
39 | std::vector< std::vector<int> > pairs(0);
|
---|
40 | std::vector<int> testPair(0);
|
---|
41 | int maxClustersInPair = _maxPairSize;
|
---|
42 | if(!maxClustersInPair)
|
---|
43 | maxClustersInPair = stableCones.size();
|
---|
44 | addClustersToPairs(testPair,pairs,distanceOK,maxClustersInPair);
|
---|
45 |
|
---|
46 | // Loop over all combinations. Calculate MidPoint. Make midPointClusters.
|
---|
47 | bool reduceConeSize = false;
|
---|
48 | for(unsigned int iPair = 0; iPair < pairs.size(); iPair++){
|
---|
49 | // Calculate rapidity, phi and pT of MidPoint.
|
---|
50 | LorentzVector midPoint(0,0,0,0);
|
---|
51 | for(unsigned int iPairMember = 0; iPairMember < pairs[iPair].size(); iPairMember++)
|
---|
52 | midPoint.add(stableCones[pairs[iPair][iPairMember]].fourVector);
|
---|
53 | iterateCone(midPoint.y(),midPoint.phi(),midPoint.pt(),towers,stableCones,reduceConeSize);
|
---|
54 | }
|
---|
55 |
|
---|
56 | sort(stableCones.begin(),stableCones.end(),ClusterPtGreater());
|
---|
57 | }
|
---|
58 |
|
---|
59 | void MidPointAlgorithm::iterateCone(double startRapidity, double startPhi, 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 | double endRapidity = trialCone.fourVector.y();
|
---|
87 | double endPhi = trialCone.fourVector.phi();
|
---|
88 | 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 |
|
---|
117 | void MidPointAlgorithm::addClustersToPairs(std::vector<int>& testPair, std::vector< std::vector<int> >& pairs,
|
---|
118 | std::vector< std::vector<bool> >& distanceOK, unsigned 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(unsigned 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(unsigned 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 |
|
---|
147 | void 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 | // Start with the highest pt cone.
|
---|
154 | std::vector<Cluster>::iterator stableConeIter1 = stableCones.begin();
|
---|
155 | if(stableConeIter1 == stableCones.end()) // Stable cone list empty?
|
---|
156 | mergingNotFinished = false;
|
---|
157 | else{
|
---|
158 | bool coneNotModified = true;
|
---|
159 | // Determine whether highest pt cone has an overlap with other stable cones.
|
---|
160 | std::vector<Cluster>::iterator stableConeIter2 = stableConeIter1;
|
---|
161 | stableConeIter2++; // 2nd highest pt cone.
|
---|
162 | while(coneNotModified && stableConeIter2 != stableCones.end()){
|
---|
163 | // Calculate overlap of the two cones.
|
---|
164 | Cluster overlap;
|
---|
165 | for(std::vector<PhysicsTower>::iterator towerIter1 = stableConeIter1->towerList.begin();
|
---|
166 | towerIter1 != stableConeIter1->towerList.end();
|
---|
167 | towerIter1++){
|
---|
168 | bool isInCone2 = false;
|
---|
169 | for(std::vector<PhysicsTower>::iterator towerIter2 = stableConeIter2->towerList.begin();
|
---|
170 | towerIter2 != stableConeIter2->towerList.end();
|
---|
171 | towerIter2++)
|
---|
172 | if(towerIter1->isEqual(*towerIter2))
|
---|
173 | isInCone2 = true;
|
---|
174 | if(isInCone2)
|
---|
175 | overlap.addTower(*towerIter1);
|
---|
176 | }
|
---|
177 | if(overlap.size()){ // non-empty overlap
|
---|
178 | coneNotModified = false;
|
---|
179 | if(overlap.fourVector.pt() >= _overlapThreshold*stableConeIter2->fourVector.pt()){
|
---|
180 | // Merge the two cones.
|
---|
181 | for(std::vector<PhysicsTower>::iterator towerIter2 = stableConeIter2->towerList.begin();
|
---|
182 | towerIter2 != stableConeIter2->towerList.end();
|
---|
183 | towerIter2++){
|
---|
184 | bool isInOverlap = false;
|
---|
185 | for(std::vector<PhysicsTower>::iterator overlapTowerIter = overlap.towerList.begin();
|
---|
186 | overlapTowerIter != overlap.towerList.end();
|
---|
187 | overlapTowerIter++)
|
---|
188 | if(towerIter2->isEqual(*overlapTowerIter))
|
---|
189 | isInOverlap = true;
|
---|
190 | if(!isInOverlap)
|
---|
191 | stableConeIter1->addTower(*towerIter2);
|
---|
192 | }
|
---|
193 | // Remove the second cone.
|
---|
194 | stableCones.erase(stableConeIter2);
|
---|
195 | }
|
---|
196 | else{
|
---|
197 | // Separate the two cones.
|
---|
198 | // Which particle goes where?
|
---|
199 | std::vector<PhysicsTower> removeFromCone1,removeFromCone2;
|
---|
200 | for(std::vector<PhysicsTower>::iterator towerIter = overlap.towerList.begin();
|
---|
201 | towerIter != overlap.towerList.end();
|
---|
202 | towerIter++){
|
---|
203 | double towerRapidity = towerIter->fourVector.y();
|
---|
204 | double towerPhi = towerIter->fourVector.phi();
|
---|
205 | // Calculate distance from cone 1.
|
---|
206 | double dRapidity1 = fabs(towerRapidity - stableConeIter1->fourVector.y());
|
---|
207 | double dPhi1 = fabs(towerPhi - stableConeIter1->fourVector.phi());
|
---|
208 | if(dPhi1 > M_PI)
|
---|
209 | dPhi1 = 2*M_PI - dPhi1;
|
---|
210 | double dRJet1 = sqrt(dRapidity1*dRapidity1 + dPhi1*dPhi1);
|
---|
211 | // Calculate distance from cone 2.
|
---|
212 | double dRapidity2 = fabs(towerRapidity - stableConeIter2->fourVector.y());
|
---|
213 | double dPhi2 = fabs(towerPhi - stableConeIter2->fourVector.phi());
|
---|
214 | if(dPhi2 > M_PI)
|
---|
215 | dPhi2 = 2*M_PI - dPhi2;
|
---|
216 | double dRJet2 = sqrt(dRapidity2*dRapidity2 + dPhi2*dPhi2);
|
---|
217 | if(dRJet1 < dRJet2)
|
---|
218 | // Particle is closer to cone 1. To be removed from cone 2.
|
---|
219 | removeFromCone2.push_back(*towerIter);
|
---|
220 | else
|
---|
221 | // Particle is closer to cone 2. To be removed from cone 1.
|
---|
222 | removeFromCone1.push_back(*towerIter);
|
---|
223 | }
|
---|
224 | // Remove particles in the overlap region from the cones to which they have the larger distance.
|
---|
225 | for(std::vector<PhysicsTower>::iterator towerIter = removeFromCone1.begin();
|
---|
226 | towerIter != removeFromCone1.end();
|
---|
227 | towerIter++)
|
---|
228 | stableConeIter1->removeTower(*towerIter);
|
---|
229 | for(std::vector<PhysicsTower>::iterator towerIter = removeFromCone2.begin();
|
---|
230 | towerIter != removeFromCone2.end();
|
---|
231 | towerIter++)
|
---|
232 | stableConeIter2->removeTower(*towerIter);
|
---|
233 | }
|
---|
234 | }
|
---|
235 | stableConeIter2++;
|
---|
236 | }
|
---|
237 | if(coneNotModified){
|
---|
238 | // Cone 1 has no overlap with any of the other cones and can become a jet.
|
---|
239 | jets.push_back(*stableConeIter1);
|
---|
240 | stableCones.erase(stableConeIter1);
|
---|
241 | }
|
---|
242 | }
|
---|
243 | }
|
---|
244 |
|
---|
245 | sort(jets.begin(),jets.end(),ClusterPtGreater());
|
---|
246 | }
|
---|
247 |
|
---|
248 | void MidPointAlgorithm::run(std::vector<PhysicsTower>& towers, std::vector<Cluster>& jets)
|
---|
249 | {
|
---|
250 | std::vector<Cluster> stableCones;
|
---|
251 | findStableConesFromSeeds(towers,stableCones);
|
---|
252 | findStableConesFromMidPoints(towers,stableCones);
|
---|
253 | splitAndMerge(stableCones,jets);
|
---|
254 | }
|
---|