source: svn/trunk/Utilities/CDFCones/src/JetCluAlgorithm.cc@ 2

#include "Utilities/CDFCones/interface/JetCluAlgorithm.h"
#include "Utilities/CDFCones/interface/ClusterComparisons.h"
#include "Utilities/CDFCones/interface/Centroid.h"
#include <algorithm>
#include <cmath>

void JetCluAlgorithm::makeSeedTowers(std::vector<PhysicsTower>& towers, std::vector<Cluster>& seedTowers)
{
  for(int iEta = 4; iEta < 48; iEta++){
    bool seg24 = true;
    if(iEta >= 8 && iEta < 14 || iEta >= 38 && iEta < 44)
      seg24 = false;
    for(int iPhi = 0; iPhi < 24; iPhi++){
      Cluster seed;
      for(std::vector<PhysicsTower>::iterator towerIter = towers.begin(); towerIter != towers.end(); towerIter++)
        if(towerIter->iEta() == iEta &&
           (seg24 && towerIter->iPhi() == iPhi || !seg24 && (towerIter->iPhi() == 2*iPhi || towerIter->iPhi() == 2*iPhi + 1)))
          seed.addTower(*towerIter);
      if(seed.centroid.Et > _seedThreshold)
        seedTowers.push_back(seed);
    }
  }
  sort(seedTowers.begin(),seedTowers.end(),ClusterCentroidEtGreater());
}

void JetCluAlgorithm::buildPreClusters(std::vector<Cluster>& seedTowers, std::vector<PhysicsTower>& towers,
                                       std::vector<Cluster>& preClusters)
{
  std::vector<Centroid> leadingSeedTowers;
  for(std::vector<Cluster>::iterator seedTowerIter = seedTowers.begin(); seedTowerIter != seedTowers.end(); seedTowerIter++){
    bool seedTowerAddedToPreCluster = false;
    std::vector<Cluster>::iterator preClusterIter = preClusters.begin();
    std::vector<Centroid>::iterator leadingSeedTowerIter = leadingSeedTowers.begin();
    while(preClusterIter != preClusters.end() && !seedTowerAddedToPreCluster){
      double dEta = fabs(seedTowerIter->centroid.eta - leadingSeedTowerIter->eta);
      double dPhi = fabs(seedTowerIter->centroid.phi - leadingSeedTowerIter->phi);
      if(dPhi > M_PI)
        dPhi = 2*M_PI - dPhi;
      if(dEta <= _coneRadius && dPhi <= _coneRadius){
        int iEtaSeedTower = seedTowerIter->towerList.begin()->iEta();
        int iPhiSeedTower = seedTowerIter->towerList.begin()->iPhi();
        if(iEtaSeedTower >= 8 && iEtaSeedTower < 14 || iEtaSeedTower >= 38 && iEtaSeedTower < 44)
          iPhiSeedTower = iPhiSeedTower/2;
        for(std::vector<PhysicsTower>::iterator preClusterTowerIter = preClusterIter->towerList.begin();
            preClusterTowerIter != preClusterIter->towerList.end() && !seedTowerAddedToPreCluster;
            preClusterTowerIter++){
          int iEtaPreClusterTower = preClusterTowerIter->iEta();
          int iPhiPreClusterTower = preClusterTowerIter->iPhi();
          if(iEtaPreClusterTower >= 8 && iEtaPreClusterTower < 14 || iEtaPreClusterTower >= 38 && iEtaPreClusterTower < 44)
            iPhiPreClusterTower = iPhiPreClusterTower/2;
          int dIEta = abs(iEtaSeedTower - iEtaPreClusterTower);
          int dIPhi = abs(iPhiSeedTower - iPhiPreClusterTower);
          if(dIPhi > 12)
            dIPhi = 24 - dIPhi;
          int adj = dIPhi*dIPhi + dIEta*dIEta;
          if(adj <= _adjacencyCut){
            for(std::vector<PhysicsTower>::iterator seedTowerTowerIter = seedTowerIter->towerList.begin();
                seedTowerTowerIter != seedTowerIter->towerList.end();
                seedTowerTowerIter++)
              preClusterIter->addTower(*seedTowerTowerIter);
            seedTowerAddedToPreCluster = true;
          }
        }
      }
      preClusterIter++;
      leadingSeedTowerIter++;
    }
    if(!seedTowerAddedToPreCluster){
      Cluster newPreCluster;
      for(std::vector<PhysicsTower>::iterator seedTowerTowerIter = seedTowerIter->towerList.begin();
          seedTowerTowerIter != seedTowerIter->towerList.end();
          seedTowerTowerIter++)
        newPreCluster.addTower(*seedTowerTowerIter);
      preClusters.push_back(newPreCluster);
      leadingSeedTowers.push_back(Centroid(newPreCluster.centroid.Et,newPreCluster.centroid.eta,newPreCluster.centroid.phi));
    }
  }
}

void JetCluAlgorithm::findStableCones(std::vector<Cluster>& preClusters, std::vector<PhysicsTower>& towers,
                                      std::vector<Cluster>& stableCones)
{
  for(std::vector<Cluster>::iterator preClusterIter = preClusters.begin(); preClusterIter != preClusters.end(); preClusterIter++){
    double startEt  = preClusterIter->centroid.Et;
    double startEta = preClusterIter->centroid.eta;
    double startPhi = preClusterIter->centroid.phi;
    int nIterations = 0;
    Cluster trialCone;
    while(nIterations++ < _maxIterations){
      trialCone.clear();
      for(std::vector<PhysicsTower>::iterator towerIter = towers.begin(); towerIter != towers.end(); towerIter++){
        double dEta = fabs(towerIter->eta() - startEta);
        double dPhi = fabs(towerIter->phi() - startPhi);
        if(dPhi > M_PI)
          dPhi = 2*M_PI - dPhi;
        double dR = sqrt(dEta*dEta + dPhi*dPhi);
        if(dR < _coneRadius)
          trialCone.addTower(*towerIter);
      }
      if(_iratch != 0)
        for(std::vector<PhysicsTower>::iterator preClusterTowerIter = preClusterIter->towerList.begin();
            preClusterTowerIter != preClusterIter->towerList.end();
            preClusterTowerIter++){
          bool foundInTrialCone = false;
          for(std::vector<PhysicsTower>::iterator trialConeTowerIter = trialCone.towerList.begin();
              trialConeTowerIter != trialCone.towerList.end() && !foundInTrialCone;
              trialConeTowerIter++)
            if(trialConeTowerIter->isEqual(*preClusterTowerIter))
              foundInTrialCone = true;
          if(!foundInTrialCone)
            trialCone.addTower(*preClusterTowerIter);
        }
      if(nIterations <= _maxIterations){
        double endEt  = trialCone.centroid.Et;
        double endEta = trialCone.centroid.eta;
        double endPhi = trialCone.centroid.phi;
        if(endEt == startEt && endEta == startEta && endPhi == startPhi)
          nIterations = _maxIterations;
        else{
          startEt  = endEt;
          startEta = endEta;
          startPhi = endPhi;
        }
      }
    }
//    bool foundIdentical = false;
//    for(std::vector<Cluster>::iterator stableConeIter = stableCones.begin();
//      stableConeIter != stableCones.end() && !foundIdentical;
//      stableConeIter++)
//      if(trialCone.centroid.isEqual(stableConeIter->centroid))
//      foundIdentical = true;
//    if(!foundIdentical)
      stableCones.push_back(trialCone);
  }
  sort(stableCones.begin(),stableCones.end(),ClusterCentroidEtGreater());
}

void JetCluAlgorithm::splitAndMerge(std::vector<Cluster>& stableCones, std::vector<Cluster>& jets)
{
  std::vector<bool> isActive;
  for(std::vector<Cluster>::iterator stableConeIter = stableCones.begin(); stableConeIter != stableCones.end(); stableConeIter++)
    isActive.push_back(bool(true));
  std::vector<bool>::iterator isActiveIter1 = isActive.begin();
  for(std::vector<Cluster>::iterator stableConeIter1 = stableCones.begin();
      stableConeIter1 != stableCones.end();
      stableConeIter1++, isActiveIter1++){
    std::vector<Cluster>::iterator stableConeIter2 = stableCones.begin();
    std::vector<bool>::iterator isActiveIter2 = isActive.begin();
    while(stableConeIter2 != stableConeIter1 && *isActiveIter1){
      if(*isActiveIter2){
        Cluster overlap;
        for(std::vector<PhysicsTower>::iterator towerIter1 = stableConeIter1->towerList.begin();
            towerIter1 != stableConeIter1->towerList.end();
            towerIter1++)
          for(std::vector<PhysicsTower>::iterator towerIter2 = stableConeIter2->towerList.begin();
              towerIter2 != stableConeIter2->towerList.end();
              towerIter2++)
            if(towerIter1->isEqual(*towerIter2)){
              overlap.addTower(*towerIter1);
              break;
            }
        if(overlap.size()){
          if(overlap.size() == stableConeIter2->size())
            *isActiveIter2 = false;
          else if(overlap.size() == stableConeIter1->size())
            *isActiveIter1 = false;
          else if(overlap.centroid.Et > _overlapThreshold*stableConeIter1->centroid.Et ||
                  overlap.centroid.Et > _overlapThreshold*stableConeIter2->centroid.Et){
            for(std::vector<PhysicsTower>::iterator stableConeTowerIter2 = stableConeIter2->towerList.begin();
                stableConeTowerIter2 != stableConeIter2->towerList.end();
                stableConeTowerIter2++){
              bool isInOverlap = false;
              for(std::vector<PhysicsTower>::iterator overlapTowerIter = overlap.towerList.begin();
                  overlapTowerIter != overlap.towerList.end() && !isInOverlap;
                  overlapTowerIter++)
                if(stableConeTowerIter2->isEqual(*overlapTowerIter))
                  isInOverlap = true;
              if(!isInOverlap)
                stableConeIter1->addTower(*stableConeTowerIter2);
            }
            *isActiveIter2 = false;
          }
          else{
            Cluster removeFromStableCone1,removeFromStableCone2,oldRemoveFromStableCone1,oldRemoveFromStableCone2;
            double etaStableCone1 = stableConeIter1->centroid.eta;
            double phiStableCone1 = stableConeIter1->centroid.phi;
            double etaStableCone2 = stableConeIter2->centroid.eta;
            double phiStableCone2 = stableConeIter2->centroid.phi;
            double dRstableCone1,dRstableCone2;
            int iterCount = 0;
            while(iterCount++ <= _maxIterations){
              oldRemoveFromStableCone1.clear();
              oldRemoveFromStableCone2.clear();
              if(iterCount > 1){
                if(removeFromStableCone1.size()){
                  Centroid stableConeCentroid1(stableConeIter1->centroid);
                  Centroid removeCentroid1(removeFromStableCone1.centroid);
                  stableConeCentroid1.subtract(removeCentroid1);
                  etaStableCone1 = stableConeCentroid1.eta;
                  phiStableCone1 = stableConeCentroid1.phi;
                }
                else{
                  etaStableCone1 = stableConeIter1->centroid.eta;
                  phiStableCone1 = stableConeIter1->centroid.phi;
                }
                if(removeFromStableCone2.size()){
                  Centroid stableConeCentroid2(stableConeIter2->centroid);
                  Centroid removeCentroid2(removeFromStableCone2.centroid);
                  stableConeCentroid2.subtract(removeCentroid2);
                  etaStableCone2 = stableConeCentroid2.eta;
                  phiStableCone2 = stableConeCentroid2.phi;
                }
                else{
                  etaStableCone2 = stableConeIter2->centroid.eta;
                  phiStableCone2 = stableConeIter2->centroid.phi;
                }
                for(std::vector<PhysicsTower>::iterator removeTowerIter1 = removeFromStableCone1.towerList.begin();
                    removeTowerIter1 != removeFromStableCone1.towerList.end();
                    removeTowerIter1++)
                  oldRemoveFromStableCone1.addTower(*removeTowerIter1);
                for(std::vector<PhysicsTower>::iterator removeTowerIter2 = removeFromStableCone2.towerList.begin();
                    removeTowerIter2 != removeFromStableCone2.towerList.end();
                    removeTowerIter2++)
                  oldRemoveFromStableCone2.addTower(*removeTowerIter2);
              }
              removeFromStableCone1.clear();
              removeFromStableCone2.clear();
              for(std::vector<PhysicsTower>::iterator overlapTowerIter = overlap.towerList.begin();
                  overlapTowerIter != overlap.towerList.end();
                  overlapTowerIter++){
                double dEta1 = fabs(overlapTowerIter->eta() - etaStableCone1);
                double dPhi1 = fabs(overlapTowerIter->phi() - phiStableCone1);
                if(dPhi1 > M_PI)
                  dPhi1 = 2*M_PI - dPhi1;
                dRstableCone1 = dEta1*dEta1 + dPhi1*dPhi1;
                double dEta2 = fabs(overlapTowerIter->eta() - etaStableCone2);
                double dPhi2 = fabs(overlapTowerIter->phi() - phiStableCone2);
                if(dPhi2 > M_PI)
                  dPhi2 = 2*M_PI - dPhi2;
                dRstableCone2 = dEta2*dEta2 + dPhi2*dPhi2;
                if(dRstableCone1 < dRstableCone2)
                  removeFromStableCone2.addTower(*overlapTowerIter);
                else
                  removeFromStableCone1.addTower(*overlapTowerIter);
              }
              if(iterCount > 1 &&
                 removeFromStableCone1.size() == oldRemoveFromStableCone1.size() &&
                 removeFromStableCone2.size() == oldRemoveFromStableCone2.size() &&
                 (!removeFromStableCone1.size() || !removeFromStableCone2.size() ||
                  removeFromStableCone1.centroid.isEqual(oldRemoveFromStableCone1.centroid) &&
                  removeFromStableCone2.centroid.isEqual(oldRemoveFromStableCone2.centroid)))
                iterCount = _maxIterations + 1;
            }
            for(std::vector<PhysicsTower>::iterator removeTowerIter1 = removeFromStableCone1.towerList.begin();
                removeTowerIter1 != removeFromStableCone1.towerList.end();
                removeTowerIter1++)
              stableConeIter1->removeTower(*removeTowerIter1);
            for(std::vector<PhysicsTower>::iterator removeTowerIter2 = removeFromStableCone2.towerList.begin();
                removeTowerIter2 != removeFromStableCone2.towerList.end();
                removeTowerIter2++)
              stableConeIter2->removeTower(*removeTowerIter2);
          }
          overlap.clear();
        }
      }
      stableConeIter2++;
      isActiveIter2++;
    }
  }
  jets.clear();
  std::vector<bool>::iterator isActiveIter = isActive.begin();
  for(std::vector<Cluster>::iterator stableConeIter = stableCones.begin();
      stableConeIter != stableCones.end();
      stableConeIter++, isActiveIter++)
    if(*isActiveIter)
      jets.push_back(*stableConeIter);
  sort(jets.begin(),jets.end(),ClusterFourVectorEtGreater());
}

void JetCluAlgorithm::run(std::vector<PhysicsTower>& towers, std::vector<Cluster>& jets)
{
  std::vector<Cluster> seedTowers;
  makeSeedTowers(towers,seedTowers);
  std::vector<Cluster> preClusters;
  buildPreClusters(seedTowers,towers,preClusters);
  std::vector<Cluster> stableCones;
  findStableCones(preClusters,towers,stableCones);
  splitAndMerge(stableCones,jets);
}