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 "AtlfastUtils/IsAssociated.h"
#include "AtlfastAlgs/GlobalEventData.h"

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

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

#include <cmath> 
#include <algorithm>
#include <assert.h>

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

#include "CLHEP/Units/SystemOfUnits.h"


namespace Atlfast {
  using std::partial_sort;


  //--------------------------------
  // Constructors and destructors
  //--------------------------------
  
  Isolator::Isolator 
  ( const std::string& name, ISvcLocator* pSvcLocator ) 
    : Algorithm( name, pSvcLocator )
    , m_tesIO( 0 ) {
    
    // Setting the parameter defaults.
    m_rClusterMatch      = 0.150;
    m_rClusterIsolation  = 0.400;
    m_eClusterIsolation  = 0.0*GeV;
    m_rCellIsolation     = 0.200;
    m_eCellIsolation     = 10.0*GeV;

    m_rLowerHalo         = 0.2;
    m_rHigherHalo        = 0.4;
    
    // Default paths for entities in the TES
    m_inputLocation                 = "/Event/AtlfastElectrons";
    m_cellLocation                  = "/Event/AtlfastCells";
    m_clusterLocation               = "/Event/AtlfastClusters";
    m_isolatedOutputLocation        = "/Event/AtlfastIsolatedElectrons";
    m_nonIsolatedOutputLocation     = "/Event/AtlfastNonIsolatedElectrons";
    
    // 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 ) ;
    
    declareProperty( "RLowerHalo",                m_rLowerHalo ) ;
    declareProperty( "RHigherHalo",               m_rHigherHalo ) ;

  }
  
  // Destructor
  Isolator::~Isolator() {
    MsgStream log( messageService(), name() ) ;
    log << MSG::INFO << "destructor" << endreq;
    delete m_tesIO;
  } 
  
  
  //---------------------------------
  // initialise() 
  //---------------------------------
  
  StatusCode Isolator::initialize(){
    MsgStream log( messageService(), name() ) ;
    log << MSG::DEBUG << "Initialising" << endreq;
    
    
    //get the Global Event Data using singleton pattern
    GlobalEventData* ged = GlobalEventData::Instance();
    m_visibleToCal   = ged -> visibleToCal();
    m_mcLocation     = ged -> mcLocation();

    //    m_tesIO = new TesIO(eventDataService());
    m_tesIO = new TesIO(m_mcLocation, ged->justHardScatter());

    
    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_nonIsolatedOutputLocation);
    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;
    const ReconstructedParticleCollection* particles(0);
    if(!m_tesIO->getDH(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<ITwoCptCell*> cells;
    if(!m_tesIO->copy<ITwoCptCellCollection> (cells, m_cellLocation)){
      
      log << MSG::DEBUG 
          << "No Cells found in TES at " 
          << m_cellLocation 
          << endreq ;
    }
    
    
    std::vector<ICluster*> clusters;
    if(!m_tesIO->copy<IClusterCollection>(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::const_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 tis;
    tis = m_tesIO->store(isolatedParticles, m_isolatedOutputLocation );  
    if( tis.isNotValid() ) {
      log << MSG::ERROR << "Failed to register output in TES at " 
          << m_isolatedOutputLocation <<  endreq ;
      return tis;
    } else{
      log << MSG::DEBUG << "Written " << isolatedParticles->size() 
          << " isolated particles to TES " << endreq;
    }
    
    tis = m_tesIO->store(nonIsolatedParticles, m_nonIsolatedOutputLocation );
    if( tis.isNotValid() ) {
      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 tis ;
    
  }
  
  
  //-------------------------------------------
  // private: isIsolated
  //------------------------------------------
  
  // Checks whether a candidate electron is isolated from clusters
  // and cells
  
  bool Isolator::isIsolated
  ( MsgStream& log, 
    ReconstructedParticleCollection::const_iterator particle,
    std::vector<ITwoCptCell*>& storedCells,
    std::vector<ICluster*>& 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 ;
    ICluster* associatedCluster = NULL ;
    
    
    
    //....................................
    // Cluster isolation
    
    if( !storedClusters.empty() ) {
      
      // Make a local copy as we want to mutate the list of pointers
      std::vector<ICluster*>  clusters( 
                                      storedClusters.begin(), 
                                      storedClusters.end() 
                                      );
      
      //Only use clusters not already claimed by ReconstructedParticles
      std::vector<ICluster*>::iterator firstNonAssociatedCluster ;
      firstNonAssociatedCluster = 
        std::partition(
                       clusters.begin(),
                       clusters.end(),
                       IsAssociated<ReconstructedParticle>
                       );
      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.empty() ) {
      
      // 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 " << endreq ;
      
      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 ) {
      
      // Calculate halo quantities for Particle
      // Scalar sum
      (*particle)->set_halo(m_kinehelp.sumETInHalo(storedCells.begin(),
                                                   storedCells.end(),
                                                   *particle,
                                                   m_rLowerHalo,
                                                   m_rHigherHalo
                                                   ));
      // Vector sum
      (*particle)->set_halovector(m_kinehelp.sumMomentumInHalo(storedCells.begin(),
                                                               storedCells.end(),
                                                               *particle,
                                                               m_rLowerHalo,
                                                               m_rHigherHalo
                                                               ));
      
      // Set the association from Particle -> Cluster
      IAOO* iamp = *particle;
      iamp->associate( associatedCluster ) ;
      
      log  << MSG::DEBUG 
           << "Particle of type "<< (*particle)->pdg_id()
           <<" associated to Cluster" 
           << endreq;
      //assert(IsAssociated<Cluster>( *particle ));
      
      // Set the association from Cluster -> Particle
      IAOO* iamc = associatedCluster;
      iamc->associate( *particle ) ;
      
      log  << MSG::DEBUG 
           << "Cluster associated to Particle of type "
           << (*particle)->pdg_id() 
           << endreq;
      
      //assert(IsAssociated<ReconstructedParticle>( associatedCluster ));
    }
    
    log << endreq ;
    
    return true ;
  }
  
} // end of namespace bracket

---