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