ClusterConeStrategy.cxx

Go to the documentation of this file.

// ================================================
//        ClusterConeStrategy class Implementation
//
//  + PreCluster helper class
//
// ================================================
//
// Namespace Atlfast
//
#include "AtlfastAlgs/ClusterConeStrategy.h"

#include "AtlfastEvent/IKinematic.h"
#include "AtlfastEvent/Cluster.h"
#include "AtlfastEvent/CollectionDefs.h"

#include "AtlfastUtils/FunctionObjects.h"

#include <cmath> 
#include <algorithm>
#include <vector>

// Gaudi includes
#include "GaudiKernel/MsgStream.h"

//Other
#include "CLHEP/Vector/LorentzVector.h"


namespace Atlfast {
  using std::sort;
  using std::partition;

  void 
  ClusterConeStrategy::makeClusters( 
                                    MsgStream& log, 
                                    const std::vector<IKinematic*>&storedCells,
                                    IKinematicVector& unusedCells,
                                    IClusterCollection* clusters) const {
    log << MSG::DEBUG << storedCells.size() << "cells to start with"<<endreq;
    // The cell collection suppplied is not mutable as it contains items
    // (or pointers to items) in the TES. 
    // Therefore we copy pointers to two local mutable collections  as we want 
    // to manipulate them as part of the clustering strategy
    
    // This is to hold cells remaining available for cluster formation. 
    localCellCollection availableCells(storedCells.begin(), 
                                       storedCells.end() );
    localCellIterator firstAvailableCell = availableCells.begin() ;
    localCellIterator lastAvailableCell  = availableCells.end() ;
    
    // This is to keep remaining candidate initiators.
    localCellCollection availableInitiators( storedCells.begin(), 
                                             storedCells.end() ) ;
    localCellIterator firstAvailableInitiator  = availableInitiators.begin() ;
    localCellIterator lastAvailableInitiator   = availableInitiators.end() ;

    // Partition the inititator candidates to retain only those above threshold
    lastAvailableInitiator 
      = partition( firstAvailableInitiator, 
                   lastAvailableInitiator,
                   PartitionCondition::AboveThresholdET( m_minInitiatorET )
                   );
    
    // Sort the remaining initiator candidates into descending order of eT
    // This means the first available initiator will always be the highest eT
    sort( firstAvailableInitiator, 
          lastAvailableInitiator, 
          SortAttribute::DescendingET() 
          );
    
    
    // Iterate over all initiator candidates
    
    while( firstAvailableInitiator != lastAvailableInitiator  ) {
      
      // Partition the cells to find those within the required R-cone
      
      //substitute an rCone with no args to stop the compiler wingeing
      // about unused parameters
      // double rCone = this->rCone( *firstAvailableInitiator ) ;
      double rCone = this->rCone( ) ;
      
      localCellIterator endCellInCone = partition
        ( 
         firstAvailableCell, 
         lastAvailableCell, 
         PartitionCondition::BelowThresholdDeltaR( *firstAvailableInitiator, 
                                                   rCone )
         );
      
      // Accumulate all cells in cone to get weighted kinematic centre
      // We use a PreCluster helper object for this
      PreCluster preclus( firstAvailableCell, endCellInCone, m_masslessJets);
      
      // Test if sum eT is above minumum. 
      if( preclus.eT() > m_minClusterET ){
        // Make a cluster and add to the collection
        HepLorentzVector temp=preclus;
        ICluster* newclus =  
          new Cluster( temp, firstAvailableCell, endCellInCone ) ;
        //      Cluster* newclus =  
        //        new Cluster( preclus, firstAvailableCell, endCellInCone ) ;
        clusters->push_back( newclus ) ;
        
        //log << MSG::DEBUG << "New Cluster created " << newclus << endreq ;
        
        // Remove all used cells from the available initiator range
        localCellIterator cell= firstAvailableCell; 
        for(; cell != endCellInCone; ++cell ){ 
          lastAvailableInitiator =
            std::remove( firstAvailableInitiator, 
                         lastAvailableInitiator, *cell );
        }
        
        // Remove all used cells from the available cell range
        firstAvailableCell = endCellInCone ;
        
      }else{
        //log << MSG::DEBUG << "Failed to form new cluster" 
        //    << *firstAvailableInitiator << endreq ;
        
        // This initiator didnt lead to a good cluster.
        // Remove it from the available initiator  range
        ++firstAvailableInitiator ;
      }
      
    }

    localCellIterator cellIter = firstAvailableCell;
    for(; cellIter!=lastAvailableCell; ++cellIter){
      //      unusedCells.push_back(new Cell(**cellIter));
      //PS 24/01/03
      //the next line replaces the previous line
      // but why was it there in the first place?
      //      unusedCells.push_back( (*cellIter)->clone() );
      unusedCells.push_back(*cellIter);
    }
    return;
  }
  
  //-------------------------------------------// private: rcone
  //------------------------------------------
  
  // looks up R-cone size according to position of initiator
  
  //double ClusterConeStrategy::rCone( Cell* initiator )
  double ClusterConeStrategy::rCone() const {
    // Temporary: return barrel always as we dont have mechanism
    // to determine calorimeter geometry yet.
    return m_rConeBarrel ;
  }
  
  
  
  
  //-------------------------------------
  // PreClusterHelper class
  //-------------------------------------
  
  // Constructor
  ClusterConeStrategy::PreCluster::PreCluster( 
                                              localCellIterator start, 
                                              localCellIterator end,
                                              bool masslessJets) : 
    m_eT_total(0),
    m_eta_weighted(0),
    m_masslessJets(masslessJets),
    m_momentum_sum( 0, 0, 0, 0 ){ 
    
    for( localCellIterator cell = start; cell != end; ++cell ){
      // Construct sum kinematic quantities weighted by eT
      // These are done EXACTLY as per the old ATLFAST++
      m_eT_total     += (*cell)->eT() ;
      m_eta_weighted += (*cell)->eta() * (*cell)->eT() ;
      
      // This next code for phi DOES NOT do things in exactly the same
      // way as old ATLFAST++
      // In particular the simple weighted phi which was implemented
      // would seem to have a problem at the cyclic boundary 
      // (unless I misunderstood, if it averaged pi-epsilon and -pi+epsilon
      // it would get 0 instead of pi
      // What is implemented here is a more generic solution which
      // is currently only used for phi, but could be used for eT
      // and eta also later
      
      m_momentum_sum += (*cell)->momentum() ;
      
    } 
  }       
  
  // queries
  
  double ClusterConeStrategy::PreCluster::eT()  { return m_eT_total ; }
  
  double ClusterConeStrategy::PreCluster::eta() { 
    // Implemented exactly as per the old ATLFAST++
    if( m_eT_total > 0 ) return m_eta_weighted / m_eT_total ;
    else return 0. ;
  }
  
  double ClusterConeStrategy::PreCluster::phi() { 
    // Implemented differently to the old ATLFAST++
    return m_momentum_sum.phi() ; 
  }
  
  
  // Conversion of eT, eta, phi to HepLorentzVector
  ClusterConeStrategy::PreCluster::operator HepLorentzVector () {
    
    // Note:  there is a good reason this doesnt use the member
    // 4-vector. It is trying to mimic as far as possible the old
    // ATLFAST++. 
    // This may be changed in the future.
    
    double theta = atan( exp( - this->eta() ) ) * 2.0 ;
    
    double x=0.,y=0.,z=0.,t=0. ;
    
    if( (0. < theta) && ( theta < M_PI ) )
      {           
        if(m_masslessJets){
          t = this->eT() / sin( theta ) ;
          x = this->eT() * cos( this->phi() ) ;
          y = this->eT() * sin( this->phi() ) ;
          z = t * cos( theta ) ;
        }else{
          t = m_momentum_sum.e();
          x = m_momentum_sum.px();
          y = m_momentum_sum.py();
          z = m_momentum_sum.pz();
        }
      }
    
    return HepLorentzVector( x, y, z, t ) ; 
  }
} // end of namespace bracket

---