---

SharedConeStrategy.cxx

Go to the documentation of this file.

#include "AtlfastAlgs/SharedConeStrategy.h"
#include "AtlfastAlgs/ClusterConeStrategy.h"
#include "AtlfastAlgs/AssocTypeRecoverer.h"
#include "AtlfastEvent/Cell.h"
#include "AtlfastEvent/Cluster.h"
#include "AtlfastEvent/CollectionDefs.h"

#include "AtlfastUtils/FunctionObjects.h"

// Gaudi includes
#include "GaudiKernel/DataSvc.h"
#include "StoreGate/DataHandle.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"

#include <assert.h>
#include <algorithm>
#include <numeric>
#include <iostream>
namespace Atlfast{

  SharedConeStrategy::SharedConeStrategy(double rConeBarrel,
                                         double rConeForward,
                                         double minInitiatorET, 
                                         double minClusterET,
                                         double barrelForwardEta): 
    m_rConeBarrel(rConeBarrel),
    m_rConeForward(rConeForward),
    m_minClusterET(minClusterET),
    m_barrelForwardEta(barrelForwardEta){

    m_coneStrategy = new ClusterConeStrategy( rConeBarrel,   
                                              rConeForward,
                                              minInitiatorET,
                                              minInitiatorET,//crazy!!!!!
                                              true);
  }

  SharedConeStrategy::~SharedConeStrategy(){delete m_coneStrategy;}

  //Function Object used by STL in makeClusters()
  class EsumIK{
  public:
    double operator()(double sum, IKinematic* ik){
      return sum + ( ik->momentum().e() );
    }
  };
  //Function Object used by STL in makeClusters()
  class CellAssociatedClusterEnergy{
  public:
    CellAssociatedClusterEnergy( 
                                ClusterCollection* clusters,
                                std::map<Cluster*,double>& cluMap,
                                double rConeBarrel
                                ):
      m_clusters(clusters), m_cluMap(cluMap), m_rConeBarrel(rConeBarrel){
    }
    void operator()(IKinematic* cell){
      
      
      // find all clusters within dR of the cell
      ClusterCollection::iterator endAssClu = 
        partition(m_clusters->begin(), 
                  m_clusters->end(), 
                  PartitionCondition::BelowThresholdDeltaR(*cell,m_rConeBarrel)
                  );
      
      //    assert(endAssClu != cluStart);
      
      //Find the sum of energies of all clusters within dR 
      // - will be used to calculate cluster weight
      //
      // in conversion from Fortran: use energy instead of 
      // eT*cosh(pseudoRapidity)
      
      double clusterEsum = 
        std::accumulate(m_clusters->begin(), endAssClu, 0., EsumIK() );
            
      //calculate the shared energy fractions
      
      ClusterCollection::iterator itr = m_clusters->begin();
      for(; itr != endAssClu; ++itr){
        
        double eCell         = cell->momentum().e();
        double eCluster      = ( *itr)->momentum().e();
        m_cluMap[*itr]      += eCell*eCluster/clusterEsum; 
        
      }

      return;
    }
  private:
    ClusterCollection* m_clusters;
    std::map<Cluster*,double>&   m_cluMap;
    double                       m_rConeBarrel;
  };
  void 
  SharedConeStrategy::makeClusters(MsgStream& log, 
                                   const std::vector<IKinematic*>& storedCells,
                                   std::vector<IKinematic*>& unusedCells,
                                   ClusterCollection* clusters) const {
    //Firstly Call the bog Standard ClusterConeStrategy
    m_coneStrategy->makeClusters(log,  storedCells, unusedCells, clusters);
    //Now recalibrate those clusters

    std::map<Cluster*,double> aMap;
    ClusterCollection::iterator cluIter = clusters->begin();
    for(; cluIter != clusters->end(); ++cluIter) {aMap[*cluIter] = 0.0;}
    
    //Make A collection of all used cells
    AssocTypeRecoverer<ClusterCollection>   atr(clusters);
    std::vector<IKinematic*>* usedCells   = atr.vectorOfIKinematics();

    // find a weight for each cluster, and put into map
    std::for_each(usedCells->begin(), 
                  usedCells->end(), 
                  CellAssociatedClusterEnergy(                 //defined above
                                              clusters,
                                              aMap,
                                              m_rConeBarrel
                                              )
                  
                  );
    delete usedCells;                       //deletes container, not constents

    //recalculate cluster energies
    cluIter = clusters->begin();
    for(; cluIter != clusters->end(); ++cluIter){
      if(aMap.find(*cluIter) != aMap.end() && aMap[*cluIter] != 0.0) {
        double theta = ( (*cluIter)->momentum() ).theta();
        double t     = aMap[*cluIter];
        double eT    = t*sin(theta); 
        
        double x = eT * cos( (*cluIter)->phi() ) ;
        double y = eT * sin( (*cluIter)->phi() ) ;
        double z =  t * cos( theta ) ;
          
        (*cluIter)->setMomentum(HepLorentzVector(x,y,z,t));
      }
    }
    //Partition away clusters below eT threshold
    ClusterCollection::iterator part = 
      std::partition(
                     clusters->begin(),
                     clusters->end(),
                     PartitionCondition::AboveThresholdET( m_minClusterET )
                     );
    //retrieve cells from discarded clusters and delete clusters
    //    std::for_each(part,clusters->end(),GetBackCells(unusedCells));
    ClusterCollection unusedClusters(part, clusters->end());
    AssocTypeRecoverer<ClusterCollection> atr2(unusedClusters);
    IKinematicCollection* recoveredCells= atr2.allAsIKinematics();
    unusedCells.reserve(unusedCells.size()+recoveredCells->size());

    IKinematicIterator iter;
    for(iter = recoveredCells->begin(); iter!= recoveredCells->end(); ++iter){
      unusedCells.push_back(*iter);
    }

    
    clusters->erase(part,clusters->end());
    //sort clusters by eT
    sort( clusters->begin(),clusters->end(),SortAttribute::DescendingET());
    return;
  }
}//end namespace

---