IslandClusterAlgo.cc

Go to the documentation of this file.

#include "RecoEcal/EgammaClusterAlgos/interface/IslandClusterAlgo.h"

// Geometry
#include "Geometry/Records/interface/IdealGeometryRecord.h"
#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
#include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
#include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"

#include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h"
#include "RecoEcal/EgammaCoreTools/interface/ClusterEtLess.h"

//
#include "DataFormats/CaloRecHit/interface/CaloID.h"

// Return a vector of clusters from a collection of EcalRecHits:
std::vector<reco::BasicCluster> IslandClusterAlgo::makeClusters(
                                  const EcalRecHitCollection* hits,
                  const CaloSubdetectorGeometry *geometry_p,
                  const CaloSubdetectorTopology *topology_p,
                  const CaloSubdetectorGeometry *geometryES_p,
                  EcalPart ecalPart,
                  bool regional,
                  const std::vector<RectangularEtaPhiRegion>& regions)
{
  seeds.clear();
  used_s.clear();
  clusters_v.clear();

  recHits_ = hits;

  double threshold = 0;
  std::string ecalPart_string;
  if (ecalPart == endcap) 
    {
      threshold = ecalEndcapSeedThreshold;
      ecalPart_string = "EndCap";
      v_chstatusSeed_ = v_chstatusSeed_Endcap_;
      v_chstatus_ = v_chstatus_Endcap_;
    }
  if (ecalPart == barrel) 
    {
      threshold = ecalBarrelSeedThreshold;
      ecalPart_string = "Barrel";
      v_chstatusSeed_ = v_chstatusSeed_Barrel_;
      v_chstatus_ = v_chstatus_Barrel_;
    }

  if (verbosity < pINFO)
    {
      std::cout << "-------------------------------------------------------------" << std::endl;
      std::cout << "Island algorithm invoked for ECAL" << ecalPart_string << std::endl;
      std::cout << "Looking for seeds, energy threshold used = " << threshold << " GeV" <<std::endl;
    }

  int nregions=0;
  if(regional) nregions=regions.size();

  if(!regional || nregions) {

    EcalRecHitCollection::const_iterator it;
    for(it = hits->begin(); it != hits->end(); it++)
      {
    double energy = it->energy();
    if (energy < threshold) continue; // need to check to see if this line is useful!

        // avoid seeding for anomalous channels
    if(! it->checkFlag(EcalRecHit::kGood)) { // if rechit is good, no need for further checks
      if (it->checkFlags( v_chstatus_ ) || it->checkFlags( v_chstatusSeed_ )) {
         continue; // the recHit has to be excluded from seeding
      }
        }

    auto thisCell = geometry_p->getGeometry(it->id());
    auto const &  position = thisCell->getPosition();

    // Require that RecHit is within clustering region in case
    // of regional reconstruction
    bool withinRegion = false;
    if (regional) {
      std::vector<RectangularEtaPhiRegion>::const_iterator region;
      for (region=regions.begin(); region!=regions.end(); region++) {
        if (region->inRegion(thisCell->etaPos(),thisCell->phiPos())) {
          withinRegion =  true;
          break;
        }
      }
    }

    if (!regional || withinRegion) {
      float ET = it->energy() * position.basicVector().unit().perp();
      if (ET > threshold) seeds.push_back(*it);
    }
      }
    
  }
  
  sort(seeds.begin(), seeds.end(), [](auto const& x, auto const& y){ return x.energy() > y.energy();});

  if (verbosity < pINFO)
    {
      std::cout << "Total number of seeds found in event = " << seeds.size() << std::endl;
    }

  mainSearch(hits,geometry_p,topology_p,geometryES_p,ecalPart);
  sort(clusters_v.rbegin(), clusters_v.rend(), isClusterEtLess);

  if (verbosity < pINFO)
    {
      std::cout << "---------- end of main search. clusters have been sorted ----" << std::endl;
    }
  
  return clusters_v; 
}


void IslandClusterAlgo::mainSearch(const EcalRecHitCollection* hits,
                                   const CaloSubdetectorGeometry *geometry_p,
                                   const CaloSubdetectorTopology *topology_p,
                                   const CaloSubdetectorGeometry *geometryES_p,
                                   EcalPart ecalPart)
{
  if (verbosity < pINFO)
    {
      std::cout << "Building clusters............" << std::endl;
    }

  // Loop over seeds:
  std::vector<EcalRecHit>::iterator it;
  for (it = seeds.begin(); it != seeds.end(); it++)
    {
      // make sure the current seed does not belong to a cluster already.
      if (used_s.find(it->id()) != used_s.end())
    {
      if (it == seeds.begin())
        {
          if (verbosity < pINFO)
        {
          std::cout << "##############################################################" << std::endl;
          std::cout << "DEBUG ALERT: Highest energy seed already belongs to a cluster!" << std::endl;
          std::cout << "##############################################################" << std::endl;
        }
        }
      continue;
    }

      // clear the vector of hits in current cluster
      current_v.clear();

      current_v.push_back( std::pair<DetId, float>(it->id(), 1.) ); // by default hit energy fractions are set at 1.
      used_s.insert(it->id());

      // Create a navigator at the seed
      CaloNavigator<DetId> navigator(it->id(), topology_p);

      searchNorth(navigator);
      navigator.home();
      searchSouth(navigator);
      navigator.home();
      searchWest(navigator, topology_p);
      navigator.home();
      searchEast(navigator, topology_p);
 
      makeCluster(hits,geometry_p,geometryES_p);
   }
}


void IslandClusterAlgo::searchNorth(const CaloNavigator<DetId> &navigator)
{
  DetId southern = navigator.pos();

  DetId northern = navigator.north();
  if (northern == DetId(0)) return; // This means that we went off the ECAL!
  // if the crystal to the north belongs to another cluster return
  if (used_s.find(northern) != used_s.end()) return;


  EcalRecHitCollection::const_iterator southern_it = recHits_->find(southern);
  EcalRecHitCollection::const_iterator northern_it = recHits_->find(northern);

  if (shouldBeAdded(northern_it, southern_it))
    {
      current_v.push_back( std::pair<DetId, float>(northern, 1.)); // by default hit energy fractions are set at 1.
      used_s.insert(northern);
      searchNorth(navigator);
    }
}


void IslandClusterAlgo::searchSouth(const CaloNavigator<DetId> &navigator)
{
  DetId northern = navigator.pos();

  DetId southern = navigator.south();
  if (southern == DetId(0)) return; // This means that we went off the ECAL!
  if (used_s.find(southern) != used_s.end()) return;


  EcalRecHitCollection::const_iterator northern_it = recHits_->find(northern);
  EcalRecHitCollection::const_iterator southern_it = recHits_->find(southern);

  if (shouldBeAdded(southern_it, northern_it))
    {
      current_v.push_back( std::pair<DetId, float>(southern, 1.)); // by default hit energy fractions are set at 1.
      used_s.insert(southern);
      searchSouth(navigator);
    }
}


void IslandClusterAlgo::searchWest(const CaloNavigator<DetId> &navigator, const CaloSubdetectorTopology* topology)
{
  DetId eastern = navigator.pos();
  EcalRecHitCollection::const_iterator eastern_it = recHits_->find(eastern);

  DetId western = navigator.west();
  if (western == DetId(0)) return; // This means that we went off the ECAL!
  EcalRecHitCollection::const_iterator western_it = recHits_->find(western);

  if (shouldBeAdded(western_it, eastern_it))
    {
      CaloNavigator<DetId> nsNavigator(western, topology);

      searchNorth(nsNavigator);
      nsNavigator.home();
      searchSouth(nsNavigator);
      nsNavigator.home();
      searchWest(navigator, topology);

      current_v.push_back( std::pair<DetId, float>(western, 1.)); // by default hit energy fractions are set at 1.
      used_s.insert(western);
    }
}


void IslandClusterAlgo::searchEast(const CaloNavigator<DetId> &navigator, const CaloSubdetectorTopology* topology)
{
  DetId western = navigator.pos();
  EcalRecHitCollection::const_iterator western_it = recHits_->find(western);

  DetId eastern = navigator.east();
  if (eastern == DetId(0)) return; // This means that we went off the ECAL!
  EcalRecHitCollection::const_iterator eastern_it = recHits_->find(eastern);

  if (shouldBeAdded(eastern_it, western_it))
    {
      CaloNavigator<DetId> nsNavigator(eastern, topology);

      searchNorth(nsNavigator);
      nsNavigator.home();
      searchSouth(nsNavigator);
      nsNavigator.home();
      searchEast(navigator, topology);

      current_v.push_back( std::pair<DetId, float>(eastern, 1.)); // by default hit energy fractions are set at 1.
      used_s.insert(eastern);
    }
}


// returns true if the candidate crystal fulfills the requirements to be added to the cluster:
bool IslandClusterAlgo::shouldBeAdded(EcalRecHitCollection::const_iterator candidate_it, EcalRecHitCollection::const_iterator previous_it)
{
  // crystal should not be included...
  if ( (candidate_it == recHits_->end())                 || // ...if it does not correspond to a hit
       (used_s.find(candidate_it->id()) != used_s.end()) || // ...if it already belongs to a cluster
       (candidate_it->energy() <= 0)                     || // ...if it has a negative or zero energy
       (candidate_it->energy() > previous_it->energy())  || // ...or if the previous crystal had lower E
       (!(candidate_it->checkFlag(EcalRecHit::kGood)) && candidate_it->checkFlags( v_chstatus_ )))
    {
      return false;
    }
  return true;
}


void IslandClusterAlgo::makeCluster(const EcalRecHitCollection* hits,
                    const CaloSubdetectorGeometry *geometry,
                    const CaloSubdetectorGeometry *geometryES)
{
  double energy = 0;
  reco::CaloID caloID;

  Point position;
  position = posCalculator_.Calculate_Location(current_v,hits,geometry,geometryES);
  
  std::vector< std::pair<DetId, float> >::iterator it;
  for (it = current_v.begin(); it != current_v.end(); it++)
    {
      EcalRecHitCollection::const_iterator itt = hits->find( (*it).first );
      EcalRecHit hit_p = *itt;
      if ( (*it).first.subdetId() == EcalBarrel ) {
              caloID = reco::CaloID::DET_ECAL_BARREL;
      } else {
              caloID = reco::CaloID::DET_ECAL_ENDCAP;
      }
      //      if (hit_p != 0)
      //    {
      energy += hit_p.energy();
      //    }
      //      else 
      //    {
      //      std::cout << "DEBUG ALERT: Requested rechit has gone missing from rechits map! :-S" << std::endl;
      //    }
    }

  if (verbosity < pINFO)
    { 
      std::cout << "******** NEW CLUSTER ********" << std::endl;
      std::cout << "No. of crystals = " << current_v.size() << std::endl;
      std::cout << "     Energy     = " << energy << std::endl;
      std::cout << "     Phi        = " << position.phi() << std::endl;
      std::cout << "     Eta        = " << position.eta() << std::endl;
      std::cout << "*****************************" << std::endl;
    }
  clusters_v.push_back(reco::BasicCluster(energy, position, caloID, current_v, reco::CaloCluster::island));
}