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