// ================================================
 //        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;
 
 /** @brief Implements the cone jet finder as a clustering strategy.
    *
    * Default R-cone value is 0.401. The 0.001 avoids problems with the
    * 0.1x0.1 grid when running on Atlfast calorimeter cells.
    * This strategy has been shown to produce different results than the JetRec
    * cone jet finder and should be superseded soon.
    */
 
   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
 
 
 
 
 
 
 
 
 
 
 

Due to the LXR bug, the updates fail sometimes to remove references to deleted files. The Saturday's full rebuilds fix these problems
This page was automatically generated by the LXR engine.