source: HiSusy/trunk/Delphes/Delphes-3.0.9/external/fastjet/plugins/CDFCones/MidPointAlgorithm.cc @ 5

Last change on this file since 5 was 5, checked in by zerwas, 11 years ago

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