SharedConeStrategy.cxx

Go to the documentation of this file.

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

#include "AtlfastUtils/FunctionObjects.h"

// Gaudi includes
#include "GaudiKernel/DataSvc.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( 
                                IClusterCollection* clusters,
                                std::map<ICluster*,double>& cluMap,
                                double rConeBarrel
                                ):
      m_clusters(clusters), m_cluMap(cluMap), m_rConeBarrel(rConeBarrel){
    }
    void operator()(const IKinematic* cell){
      
      
      // find all clusters within dR of the cell
      IClusterCollection::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
      
      IClusterCollection::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:
    IClusterCollection* m_clusters;
    std::map<ICluster*,double>&   m_cluMap;
    double                       m_rConeBarrel;
  };
  void 
  SharedConeStrategy::makeClusters(MsgStream& log, 
                                   const std::vector<IKinematic*>& storedCells,
                                   IKinematicVector& unusedCells,
                                   IClusterCollection* clusters) const {
    //Firstly Call the bog Standard ClusterConeStrategy
    m_coneStrategy->makeClusters(log,  storedCells, unusedCells, clusters);
    //Now recalibrate those clusters

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

    // find a weight for each cluster, and put into map
    std::for_each(usedCells.begin(), 
                  usedCells.end(), 
                  CellAssociatedClusterEnergy(                 //defined above
                                              clusters,
                                              aMap,
                                              m_rConeBarrel
                                              )
                  
                  );

    //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
    IClusterCollection::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));
    IClusterCollection unusedClusters(part, clusters->end());
    //    AssocTypeRecoverer atr2(unusedClusters.begin(), unusedClusters.end());
    //    IKinematicVector* recoveredCells= atr2.ikinematics();
    TypeVisitor tv = 
      ContainerAssocsDispatcher(
                               unusedClusters.begin(), 
                               unusedClusters.end(),
                               TypeVisitor()
                               );
    IKinematicVector recoveredCells = tv.ikinematics();

    unusedCells.reserve(unusedCells.size()+recoveredCells.size());

    std::copy( 
              recoveredCells.begin(), 
              recoveredCells.end(), 
              back_inserter(unusedCells) 
              );

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

---