---

Isolator.cxx

Go to the documentation of this file.

// ================================================
// Isolator class Implementation
// ================================================
//
// THIS TEXT TO BE REPLACED BY ATLAS STANDARD FORMAT
//
// Namespace ATLFast
//

#include "AtlfastAlgs/Isolator.h"
#include "AtlfastAlgs/ClusterIsAssoc.h"
#include "AtlfastAlgs/GlobalEventData.h"

#include "AtlfastEvent/Phi.h"
#include "AtlfastEvent/Cluster.h"

#include "AtlfastUtils/FunctionObjects.h"
#include "AtlfastUtils/HeaderPrinter.h"
#include "AtlfastUtils/HepMC_helper/IMCselector.h"

#include <cmath> 
#include <algorithm>

// Gaudi includes
#include "GaudiKernel/DataSvc.h"
#include "StoreGate/DataHandle.h"
#include "AtlfastEvent/MsgStreamDefs.h"

namespace Atlfast {
  using std::partial_sort;


  //--------------------------------
  // Constructors and destructors
  //--------------------------------
  
  Isolator::Isolator 
  ( const std::string& name, ISvcLocator* pSvcLocator ) 
    : Algorithm( name, pSvcLocator ){
    
    // Setting the parameter defaults.
    m_rClusterMatch      = 0.150;
    m_rClusterIsolation  = 0.400;
    m_eClusterIsolation  = 0.0;
    m_rCellIsolation     = 0.200;
    m_eCellIsolation     = 10.0;
    
    // Default paths for entities in the TES
    m_inputLocation                 = "Unknown";
    m_cellLocation                  = "/Event/AtlfastCell";
    m_clusterLocation               = "/Event/AtlfastCluster";
    m_isolatedOutputLocation        = "Unknown";
    m_nonIsolatedOutputLocation     = "Unknown";
    
    // This is how you declare the paramemters to Gaudi so that
    // they can be over-written via the job options file
    
    declareProperty( "RClusterMatch",             m_rClusterMatch ) ;
    declareProperty( "RClusterIsolation",         m_rClusterIsolation ) ;
    declareProperty( "EClusterIsolation",         m_eClusterIsolation ) ;
    declareProperty( "RCellIsolation",            m_rCellIsolation ) ;
    declareProperty( "ECellIsolation",            m_eCellIsolation ) ;
    
    declareProperty( "InputLocation",             m_inputLocation ) ;
    declareProperty( "CellLocation",              m_cellLocation ) ;
    declareProperty( "ClusterLocation",           m_clusterLocation ) ;
    declareProperty( "IsolatedOutputLocation",    m_isolatedOutputLocation ) ;
    declareProperty( "NonIsolatedOutputLocation", m_nonIsolatedOutputLocation ) ;
    
  }
  
  // Destructor
  Isolator::~Isolator() {
    MsgStream log( messageService(), name() ) ;
    log << MSG::INFO << "destructor" << endreq;
    if (m_tesIO) {
      delete m_tesIO;
    }
  } 
  
  
  //---------------------------------
  // initialise() 
  //---------------------------------
  
  StatusCode Isolator::initialize(){
    MsgStream log( messageService(), name() ) ;
    log << MSG::DEBUG << "Initialising" << endreq;
    
    
    //    m_tesIO = new TesIO(eventDataService());
    m_tesIO = new TesIO();
    //get the Global Event Data using singleton pattern
    GlobalEventData* ged = GlobalEventData::Instance();
    m_visibleToCal   = ged -> visibleToCal();
    
    HeaderPrinter hp("Atlfast Isolator:", log);
    
    hp.add("TES Locations:             ");
    hp.add(" Particles from            ", m_inputLocation); 
    hp.add(" Cells from                ", m_cellLocation);    
    hp.add(" Clusters from             ", m_clusterLocation); 
    hp.add(" Isolated Particles to     ", m_isolatedOutputLocation);
    hp.add(" Non-Isolated Particles to ", m_isolatedOutputLocation);
    hp.add("Cluster match DeltaR       ", m_rClusterMatch);
    hp.add("Cluster isolation DeltaR   ", m_rClusterIsolation);
    hp.add("Cluster isolation E thresh ", m_eClusterIsolation);
    hp.add("Cell isolation DeltaR      ", m_rCellIsolation);
    hp.add("Cell isolation E thresh    ", m_eCellIsolation);
    hp.print();
    
    
    return StatusCode::SUCCESS ;
  }
  
  //---------------------------------
  // finalise() 
  //---------------------------------
  
  StatusCode Isolator::finalize(){
    MsgStream log( messageService(), name() ) ;
    log << MSG::INFO << "finalizing" << endreq;
    return StatusCode::SUCCESS ;
  }
  
  
  //----------------------------------------------
  // execute() method called once per event
  //----------------------------------------------
  //
  //  This algorithm creates isolated electrons. 
  //  It scans the list of MC electrons. If an electron is found it creates 
  //  a candidate Electron object with a smeared energy.  Different 
  //  energy-smearing for low luminosity and high luminosity can be invoked.
  //  It then checks several isolation criteria with respect to calorimeter
  //  clusters and cells.
  //  Only if the isolation is satisfied is the Electron kept.
  //  Finally all successful ELectrons are added to the event store.
  
  
  StatusCode Isolator::execute( ){
    
    //................................
    // make a message logging stream
    
    MsgStream log( messageService(), name() ) ;
    std::string mess;
    //...............................
    // Extract the input particles which are to be tested for isolation.
    
    //std::vector<ReconstructedParticle*> particles;
    ReconstructedParticleCollection particles;
    if(!m_tesIO->
       copy<ReconstructedParticleCollection >(particles, m_inputLocation)){
      log << MSG::DEBUG 
          << "No reconstructed particles in TES to treat" 
          << endreq ;
      return StatusCode::SUCCESS ;
    }
    
    log << MSG::DEBUG 
        << "Found " 
        << particles.size() 
        << " particles in TES "
        << endreq ;
    
    
    //.........................................................
    // Extract the (pointers to) cells and clusters from the TES into 
    // locally defined arrays.
    
    
    std::vector<Cell*> cells;
    if(!m_tesIO->copy<CellCollection> (cells, m_cellLocation)){
      
      log << MSG::DEBUG 
          << "No Cells found in TES at " 
          << m_cellLocation 
          << endreq ;
    }
    
    
    std::vector<Cluster*> clusters;
    if(!m_tesIO->copy<ClusterCollection>(clusters, m_clusterLocation)){
      log << MSG::DEBUG 
          << "No Clusters found in TES at " 
          << m_clusterLocation 
          << endreq ;
    }
    
    
    // ......................
    // Make an empty container object which will store (pointers to) all 
    // sucessfully isolated and non isolated electrons.  
    
    ReconstructedParticleCollection* isolatedParticles = 
      new ReconstructedParticleCollection ;
    
    ReconstructedParticleCollection* nonIsolatedParticles = 
      new ReconstructedParticleCollection ;
    
    
    //.........................................................
    // For each particle in the input list test if it is isolated.
    // If so copy it into the output collection.
    
    ReconstructedParticleCollection::iterator particle ;
    
    for( particle = particles.begin(); 
         particle != particles.end(); 
         ++particle)
      {       
        if( this->isIsolated( log, particle, cells, clusters ) )     
          isolatedParticles->push_back
            ( new ReconstructedParticle( *particle ) ); 
        else
          nonIsolatedParticles->push_back
            ( new ReconstructedParticle( *particle ) );         
      }
    
    
    //......................................
    // Register particles in the transient event store. 
    
    TesIoStat sc;
    sc = m_tesIO->store(isolatedParticles, m_isolatedOutputLocation );  
    
    if( !sc ) {
      log << MSG::ERROR << "Failed to register output in TES at " 
          << m_isolatedOutputLocation <<  endreq ;
      return sc;
    } else{
      log << MSG::DEBUG << "Written " << isolatedParticles->size() 
          << " isolated particles to TES " << endreq;
    }
    
    sc = m_tesIO->store(nonIsolatedParticles, m_nonIsolatedOutputLocation );
    
    if( !sc ) {
      log << MSG::ERROR << "Failed to register output in TES at " 
          << m_nonIsolatedOutputLocation <<  endreq ;
    }
    else {
      log << MSG::DEBUG << "Written " << nonIsolatedParticles->size() 
          << " NON-isolated particles to TES " << endreq;
    }
    
    return sc ;
    
  }
  
  
  //-------------------------------------------
  // private: isIsolated
  //------------------------------------------
  
  // Checks whether a candidate electron is isolated from clusters
  // and cells
  
  bool Isolator::isIsolated
  ( MsgStream& log, 
    ReconstructedParticleCollection::iterator particle,
    std::vector<Cell*>& storedCells,
    std::vector<Cluster*>& storedClusters
    
    ){
    
    log << MSG::DEBUG << "\n\t Treating particle: " << *particle ;
    
    // Things we need to remember until the end of the method
    // If, during processing of Clusters, we identify an "associated cluster"
    // we remember it for use if the particle is finally judged to be isolated
    bool associated = false ;
    Cluster* associatedCluster = NULL ;
    
    
    
    //....................................
    // Cluster isolation
    
    if( storedClusters.begin() ) {
      
      if( storedClusters.size() > 0 ) {
        
        // Make a local copy as we want to mutate the list of pointers
        std::vector<Cluster*>  clusters( 
                                        storedClusters.begin(), 
                                        storedClusters.end() 
                                        );
        
        //Only use clusters not already claimed by ReconstructedParticles
        std::vector<Cluster*>::iterator firstNonAssociatedCluster ;
        firstNonAssociatedCluster = std::partition(
                                                   clusters.begin(),
                                                   clusters.end(),
                                                   ClusterIsAssoc()
                                                   );
        if(firstNonAssociatedCluster != clusters.begin()){
          log<<MSG::DEBUG<<"First "
             <<firstNonAssociatedCluster-clusters.begin()
             <<" clusters are associated to ReconstructedParticles "
             <<clusters.end()-firstNonAssociatedCluster
             <<" cluster(s) remain"
             <<endreq;
        }else{
          log<<MSG::DEBUG
             <<"No clusters are associated to ReconstructedParticles"
             <<endreq;
        }
        if(clusters.end()!=firstNonAssociatedCluster){
          // .....................
          // Associate unclaimed clusters with current particle. 
          // This looks at the nearest cluster only, and decides whether
          // it is within an "association" cone
          
          partial_sort(
                       firstNonAssociatedCluster,
                       firstNonAssociatedCluster+1,
                       clusters.end(), 
                       SortAttribute::DeltaR( *particle )
                       );
          
          log<<MSG::DEBUG<<"Sorted clusters "<<endreq;
          
          associated =
            (m_kinehelp.deltaR(*particle, *firstNonAssociatedCluster) <  
             m_rClusterMatch);
          
          if( associated ) 
            {
              log << "\n\t -->Found associated cluster: "  
                  << *clusters.begin() ;
              
              // Remember the relevant cluster until the end of this method
              associatedCluster = *firstNonAssociatedCluster ;
              
              // Set the first nonassociated cluster
              ++firstNonAssociatedCluster;
              
            }
          
          
          //............................................
          // Now we can check Cluster Isolation 
          // (if there is at least one more cluster to check)
          
          // Use kinematic helper to do all work
          double eClusterSum = m_kinehelp.sumETInCone( 
                                                      firstNonAssociatedCluster, 
                                                      clusters.end(), 
                                                      *particle,
                                                      m_rClusterIsolation
                                                      );
          
          // Apply the cut criteria
          
          if( eClusterSum  > m_eClusterIsolation ) 
            log  << "\n\t -->Excess Cluster energy in cone = " << eClusterSum 
                 << "   => NOT-ISOLATED " << endreq;       
          else
            log  << "\n\t -->Excess Cluster energy in cone = " << eClusterSum 
                 << "   => ISOLATED "  ;
          
          if( eClusterSum  > m_eClusterIsolation )  return false ;
          
        }
        
      }
    }
    
    
    //............................................
    // Cell Isolation: 
    
    if( storedCells.begin() ) {
      
      // Test on isolation in respect of the total transverese energy of all 
      // Cells within a given r-cone
      
      // Use kinematic helper to do all work
      double eCellSum = m_kinehelp.sumETInCone( 
                                               storedCells.begin(), 
                                               storedCells.end(), 
                                               *particle,
                                               m_rCellIsolation
                                               );
      
      
      // Subtract the energy associated with particle under test
      // as the cells summed will have included its deposit
      //    CalSelect depositsInCal;
      //    if(depositsInCal( *particle) ) eCellSum -= (*particle)->eT() ;
      bool itDeposits = m_visibleToCal->operator()( (*particle)->truth() );
      
      if(itDeposits) eCellSum -= (*particle)->eT() ;
      
      
      // Apply isolation criteria
      
      if( eCellSum  > m_eCellIsolation ) 
        log  
          << "\n\t -->Excess Cell energy in cone = " << eCellSum 
          << "   => NOT-ISOLATED " << endreq ;    
      else
        log   
          << "\n\t -->Excess Cell energy in cone = " << eCellSum 
          << "   => ISOLATED " ;
      
      if( eCellSum  > m_eCellIsolation ) return false ;
      
    }
    
    //..........................................
    // If we get here then the candidate was judged to be isolated
    
    // It is only now that we wish to set associations between isolated
    // particles and clusters. I.e we dont want to flag associations for 
    // non-isolated particles
    
    if( associated ) {
      
      // Set the association from Particle -> Cluster
      IAssociationManager* iamp = *particle;
      iamp->associate( associatedCluster ) ;
      
      // Set the association from Cluster -> Particle
      IAssociationManager* iamc = associatedCluster;
      iamc->associate( *particle ) ;
      
    }
    
    log << endreq ;
    
    return true ;
  }
  
} // end of namespace bracket

---