Fork me on GitHub

source: git/external/fastjet/plugins/CDFCones/MidPointAlgorithm.cc@ 45d8342

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 45d8342 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: 14.4 KB
RevLine 
[d7d2da3]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 MidPointAlgorithm.cc file
12//
[35cdc46]13// 2014-08-13 Matteo Cacciari and Gavin Salam
14// * changed a number of int -> unsigned (and in one case
15// added explicit conversion to int) to eliminate
16// long-standing compiler warnings
17//
[d7d2da3]18// 2009-01-17 Gregory Soyez <soyez@fastjet.fr>
19//
20// * put the code in the fastjet::cdf namespace
21//
22// 2008-01-15 Gavin Salam <salam@lpthe.jussieu.fr>
23//
24// * put in the volatile fix for midpoint with optimization;
25// checked regression test, and that speed comes out OK
26//
27//
28// 2007-10-16 Gavin Salam <salam@lpthe.jussieu.fr>
29//
30// * added protection to midpoint for cases where no stable cones
31// are found
32//
33// 2007-03-10 Gavin Salam <salam@lpthe.jussieu.fr>
34//
35// * added support for the pttilde scale choice in the CDF midpoint code
36//
37// 2007-02-21 Gavin Salam <salam@lpthe.jussieu.fr>
38//
39// * added option of choosing the scale used in the split-merge
40// procedure (pt [default], Et or mt)
41//
42// 2006-09-24 Gavin Salam <salam@lpthe.jussieu.fr>
43//
44// * removed the local_sort method
45//
46// 2006-09-24 Gavin Salam <salam@lpthe.jussieu.fr>
47//
48// * added JetClu+MidPoint to FastJet
49
50#include "MidPointAlgorithm.hh"
51#include "ClusterComparisons.hh"
52#include <algorithm>
53#include <iostream>
54#include <cmath>
55
56#include <fastjet/internal/base.hh>
57
58FASTJET_BEGIN_NAMESPACE
59
60namespace cdf{
61
62void MidPointAlgorithm::findStableConesFromSeeds(std::vector<PhysicsTower>& towers, std::vector<Cluster>& stableCones)
63{
64 bool reduceConeSize = true;
65 for(std::vector<PhysicsTower>::iterator towerIter = towers.begin(); towerIter != towers.end(); towerIter++)
66 if(towerIter->fourVector.pt() > _seedThreshold)
67 iterateCone(towerIter->fourVector.y(),towerIter->fourVector.phi(),0,towers,stableCones,reduceConeSize);
68}
69
70void MidPointAlgorithm::findStableConesFromMidPoints(std::vector<PhysicsTower>& towers, std::vector<Cluster>& stableCones)
71{
72 // distanceOK[i-1][j] = Is distance between stableCones i and j (i>j) less than 2*_coneRadius?
73 std::vector< std::vector<bool> > distanceOK;
74 distanceOK.resize(stableCones.size() - 1);
[35cdc46]75 // MC+GPS 2014-08-13, replaced int with unsigned
76 for(unsigned nCluster1 = 1; nCluster1 < stableCones.size(); nCluster1++){
[d7d2da3]77 distanceOK[nCluster1 - 1].resize(nCluster1);
78 double cluster1Rapidity = stableCones[nCluster1].fourVector.y();
79 double cluster1Phi = stableCones[nCluster1].fourVector.phi();
[35cdc46]80 // MC+GPS 2014-08-13, replaced int with unsigned
81 for(unsigned nCluster2 = 0; nCluster2 < nCluster1; nCluster2++){
[d7d2da3]82 double cluster2Rapidity = stableCones[nCluster2].fourVector.y();
83 double cluster2Phi = stableCones[nCluster2].fourVector.phi();
84 double dRapidity = fabs(cluster1Rapidity - cluster2Rapidity);
85 double dPhi = fabs(cluster1Phi - cluster2Phi);
86 if(dPhi > M_PI)
87 dPhi = 2*M_PI - dPhi;
88 double dR = sqrt(dRapidity*dRapidity + dPhi*dPhi);
89 distanceOK[nCluster1 - 1][nCluster2] = dR < 2*_coneRadius;
90 }
91 }
92
93 // Find all pairs (triplets, ...) of stableCones which are less than 2*_coneRadius apart from each other.
94 std::vector< std::vector<int> > pairs(0);
95 std::vector<int> testPair(0);
96 int maxClustersInPair = _maxPairSize;
97 if(!maxClustersInPair)
98 maxClustersInPair = stableCones.size();
99 addClustersToPairs(testPair,pairs,distanceOK,maxClustersInPair);
100
101 // Loop over all combinations. Calculate MidPoint. Make midPointClusters.
102 bool reduceConeSize = false;
[35cdc46]103 // MC+GPS 2014-08-13, replaced int with unsigned
104 for(unsigned iPair = 0; iPair < pairs.size(); iPair++){
[d7d2da3]105 // Calculate rapidity, phi and pT of MidPoint.
106 LorentzVector midPoint(0,0,0,0);
[35cdc46]107 // MC+GPS 2014-08-13, replaced int with unsigned
108 for(unsigned iPairMember = 0; iPairMember < pairs[iPair].size(); iPairMember++)
[d7d2da3]109 midPoint.add(stableCones[pairs[iPair][iPairMember]].fourVector);
110 iterateCone(midPoint.y(),midPoint.phi(),midPoint.pt(),towers,stableCones,reduceConeSize);
111 }
112
113 //sort(stableCones.begin(),stableCones.end(),ClusterPtGreater());
114 local_sort(stableCones); // GPS mod. to allow split-merge with various scales
115}
116
117
118void MidPointAlgorithm::iterateCone(volatile double startRapidity, volatile double startPhi, volatile double startPt,
119 std::vector<PhysicsTower>& towers, std::vector<Cluster>& stableCones, bool reduceConeSize)
120{
121 int nIterations = 0;
122 bool keepJet = true;
123 Cluster trialCone;
124 double iterationConeRadius = _coneRadius;
125 if(reduceConeSize)
126 iterationConeRadius *= sqrt(_coneAreaFraction);
127 while(nIterations++ < _maxIterations + 1 && keepJet){
128 trialCone.clear();
129 // Find particles which should go in the cone.
130 if(nIterations == _maxIterations + 1)
131 iterationConeRadius = _coneRadius;
132 for(std::vector<PhysicsTower>::iterator towerIter = towers.begin(); towerIter != towers.end(); towerIter++){
133 double dRapidity = fabs(towerIter->fourVector.y() - startRapidity);
134 double dPhi = fabs(towerIter->fourVector.phi() - startPhi);
135 if(dPhi > M_PI)
136 dPhi = 2*M_PI - dPhi;
137 double dR = sqrt(dRapidity*dRapidity + dPhi*dPhi);
138 if(dR < iterationConeRadius)
139 trialCone.addTower(*towerIter);
140 }
141 if(!trialCone.size()) // Empty cone?
142 keepJet = false;
143 else{
144 if(nIterations <= _maxIterations){
145 volatile double endRapidity = trialCone.fourVector.y();
146 volatile double endPhi = trialCone.fourVector.phi();
147 volatile double endPt = trialCone.fourVector.pt();
148 // Do we have a stable cone?
149 if(endRapidity == startRapidity && endPhi == startPhi && endPt == startPt){
150 // If cone size is reduced, then do one more iteration.
151 nIterations = _maxIterations;
152 if(!reduceConeSize)
153 nIterations++;
154 }
155 else{
156 // Another iteration.
157 startRapidity = endRapidity;
158 startPhi = endPhi;
159 startPt = endPt;
160 }
161 }
162 }
163 }
164
165 if(keepJet){
166 // We have a stable cone.
167 bool identical = false;
168 for(std::vector<Cluster>::iterator stableConeIter = stableCones.begin(); stableConeIter != stableCones.end(); stableConeIter++)
169 if(trialCone.fourVector.isEqual(stableConeIter->fourVector))
170 identical = true;
171 if(!identical)
172 stableCones.push_back(trialCone);
173 }
174}
175
176void MidPointAlgorithm::addClustersToPairs(std::vector<int>& testPair, std::vector< std::vector<int> >& pairs,
177 std::vector< std::vector<bool> >& distanceOK, int maxClustersInPair)
178{
179 // Recursively adds clusters to pairs, triplets, ... whose mid-points are then calculated.
180
181 // Find StableCone number to start with (either 0 at the beginning or last element of testPair + 1).
182 int nextClusterStart = 0;
183 if(testPair.size())
184 nextClusterStart = testPair.back() + 1;
[35cdc46]185 // MC+GPS 2014-08-13, replaced int nextCluster with unsigned
186 for(unsigned nextCluster = nextClusterStart; nextCluster <= distanceOK.size(); nextCluster++){
[d7d2da3]187 // Is new SeedCone less than 2*_coneRadius apart from all clusters in testPair?
188 bool addCluster = true;
[35cdc46]189 // MC+GPS 2014-08-13, replaced int iCluster with unsigned
190 for(unsigned iCluster = 0; iCluster < testPair.size() && addCluster; iCluster++)
[d7d2da3]191 if(!distanceOK[nextCluster - 1][testPair[iCluster]])
192 addCluster = false;
193 if(addCluster){
194 // Add it to the testPair.
195 testPair.push_back(nextCluster);
196 // If testPair is a pair, add it to pairs.
197 if(testPair.size() > 1)
198 pairs.push_back(testPair);
199 // If not bigger than allowed, find more clusters within 2*_coneRadius.
[35cdc46]200 // GPS+MC 2014-08-13, replaced testPair.size() with int(testPair.size())
201 if(int(testPair.size()) < maxClustersInPair)
[d7d2da3]202 addClustersToPairs(testPair,pairs,distanceOK,maxClustersInPair);
203 // All combinations containing testPair found. Remove last element.
204 testPair.pop_back();
205 }
206 }
207}
208
209void MidPointAlgorithm::splitAndMerge(std::vector<Cluster>& stableCones, std::vector<Cluster>& jets)
210{
211 bool mergingNotFinished = true;
212 while(mergingNotFinished){
213 // Sort the stable cones (highest pt first).
214 //sort(stableCones.begin(),stableCones.end(),ClusterPtGreater());
215 local_sort(stableCones);
216
217 // Start with the highest pt cone.
218 std::vector<Cluster>::iterator stableConeIter1 = stableCones.begin();
219 if(stableConeIter1 == stableCones.end()) // Stable cone list empty?
220 mergingNotFinished = false;
221 else{
222 bool coneNotModified = true;
223 // Determine whether highest pt cone has an overlap with other stable cones.
224 std::vector<Cluster>::iterator stableConeIter2 = stableConeIter1;
225 stableConeIter2++; // 2nd highest pt cone.
226 while(coneNotModified && stableConeIter2 != stableCones.end()){
227 // Calculate overlap of the two cones.
228 Cluster overlap;
229 for(std::vector<PhysicsTower>::iterator towerIter1 = stableConeIter1->towerList.begin();
230 towerIter1 != stableConeIter1->towerList.end();
231 towerIter1++){
232 bool isInCone2 = false;
233 for(std::vector<PhysicsTower>::iterator towerIter2 = stableConeIter2->towerList.begin();
234 towerIter2 != stableConeIter2->towerList.end();
235 towerIter2++)
236 if(towerIter1->isEqual(*towerIter2))
237 isInCone2 = true;
238 if(isInCone2)
239 overlap.addTower(*towerIter1);
240 }
241 if(overlap.size()){ // non-empty overlap
242 coneNotModified = false;
243 // GPS mod to allow various scale choices in split merge --------
244 double overlap_scale, jet2_scale;
245 switch(_smScale) {
246 case SM_pt:
247 overlap_scale = overlap.fourVector.pt();
248 jet2_scale = stableConeIter2->fourVector.pt();
249 break;
250 case SM_Et:
251 overlap_scale = overlap.fourVector.Et();
252 jet2_scale = stableConeIter2->fourVector.Et();
253 break;
254 case SM_mt:
255 overlap_scale = overlap.fourVector.mt();
256 jet2_scale = stableConeIter2->fourVector.mt();
257 break;
258 case SM_pttilde:
259 overlap_scale = overlap.pt_tilde;
260 jet2_scale = stableConeIter2->pt_tilde;
261 break;
262 default:
263 std::cerr << "Unrecognized value for _smScale: "
264 << _smScale << std::endl;
265 exit(-1);
266 }
267 if(overlap_scale >= _overlapThreshold*jet2_scale){
268 // end of GPS modification ---------------------------
269 //if(overlap.fourVector.pt() >= _overlapThreshold*stableConeIter2->fourVector.pt()){
270 // Merge the two cones.
271 for(std::vector<PhysicsTower>::iterator towerIter2 = stableConeIter2->towerList.begin();
272 towerIter2 != stableConeIter2->towerList.end();
273 towerIter2++){
274 bool isInOverlap = false;
275 for(std::vector<PhysicsTower>::iterator overlapTowerIter = overlap.towerList.begin();
276 overlapTowerIter != overlap.towerList.end();
277 overlapTowerIter++)
278 if(towerIter2->isEqual(*overlapTowerIter))
279 isInOverlap = true;
280 if(!isInOverlap)
281 stableConeIter1->addTower(*towerIter2);
282 }
283 // Remove the second cone.
284 stableCones.erase(stableConeIter2);
285 }
286 else{
287 // Separate the two cones.
288 // Which particle goes where?
289 std::vector<PhysicsTower> removeFromCone1,removeFromCone2;
290 for(std::vector<PhysicsTower>::iterator towerIter = overlap.towerList.begin();
291 towerIter != overlap.towerList.end();
292 towerIter++){
293 double towerRapidity = towerIter->fourVector.y();
294 double towerPhi = towerIter->fourVector.phi();
295 // Calculate distance from cone 1.
296 double dRapidity1 = fabs(towerRapidity - stableConeIter1->fourVector.y());
297 double dPhi1 = fabs(towerPhi - stableConeIter1->fourVector.phi());
298 if(dPhi1 > M_PI)
299 dPhi1 = 2*M_PI - dPhi1;
300 double dRJet1 = sqrt(dRapidity1*dRapidity1 + dPhi1*dPhi1);
301 // Calculate distance from cone 2.
302 double dRapidity2 = fabs(towerRapidity - stableConeIter2->fourVector.y());
303 double dPhi2 = fabs(towerPhi - stableConeIter2->fourVector.phi());
304 if(dPhi2 > M_PI)
305 dPhi2 = 2*M_PI - dPhi2;
306 double dRJet2 = sqrt(dRapidity2*dRapidity2 + dPhi2*dPhi2);
307 if(dRJet1 < dRJet2)
308 // Particle is closer to cone 1. To be removed from cone 2.
309 removeFromCone2.push_back(*towerIter);
310 else
311 // Particle is closer to cone 2. To be removed from cone 1.
312 removeFromCone1.push_back(*towerIter);
313 }
314 // Remove particles in the overlap region from the cones to which they have the larger distance.
315 for(std::vector<PhysicsTower>::iterator towerIter = removeFromCone1.begin();
316 towerIter != removeFromCone1.end();
317 towerIter++)
318 stableConeIter1->removeTower(*towerIter);
319 for(std::vector<PhysicsTower>::iterator towerIter = removeFromCone2.begin();
320 towerIter != removeFromCone2.end();
321 towerIter++)
322 stableConeIter2->removeTower(*towerIter);
323 }
324 }
325 stableConeIter2++;
326 }
327 if(coneNotModified){
328 // Cone 1 has no overlap with any of the other cones and can become a jet.
329 jets.push_back(*stableConeIter1);
330 stableCones.erase(stableConeIter1);
331 }
332 }
333 }
334
335 //sort(jets.begin(),jets.end(),ClusterPtGreater());
336 local_sort(jets); // GPS mod. to allow split-merge with various scales
337}
338
339
340
341void MidPointAlgorithm::local_sort(std::vector<Cluster>& clusters) {
342 switch(_smScale) {
343 case SM_pt:
344 sort(clusters.begin(),clusters.end(),ClusterPtGreater());
345 break;
346 case SM_Et:
347 sort(clusters.begin(),clusters.end(),ClusterFourVectorEtGreater());
348 break;
349 case SM_mt:
350 sort(clusters.begin(),clusters.end(),ClusterMtGreater());
351 break;
352 case SM_pttilde:
353 sort(clusters.begin(),clusters.end(),ClusterPtTildeGreater());
354 break;
355 default:
356 std::cerr << "Unrecognized value for _smScale: " << _smScale << std::endl;
357 exit(-1);
358 }
359}
360
361
362
363void MidPointAlgorithm::run(std::vector<PhysicsTower>& towers, std::vector<Cluster>& jets)
364{
365 std::vector<Cluster> stableCones;
366 findStableConesFromSeeds(towers,stableCones);
367 // GPS addition to prevent crashes if no stable cones
368 // are found (e.g. all particles below seed threshold)
369 if (stableCones.size() > 0) {
370 findStableConesFromMidPoints(towers,stableCones);
371 splitAndMerge(stableCones,jets);
372 }
373}
374
375} // namespace cdf
376
377FASTJET_END_NAMESPACE
Note: See TracBrowser for help on using the repository browser.