KIsolator.cxx

Go to the documentation of this file.

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

#include "AtlfastAlgs/KIsolator.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
  //--------------------------------
  
  KIsolator::KIsolator 
  ( 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;
    
    // 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 ) ;
    
  }
  
  // Destructor
  KIsolator::~KIsolator() {
    MsgStream log( messageService(), name() ) ;
    log << MSG::INFO << "destructor" << endreq;
    if (m_tesIO) delete m_tesIO;
  } 
  
  
  //---------------------------------
  // initialise() 
  //---------------------------------
  
  StatusCode KIsolator::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 KIsolator:", 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 KIsolator::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 KIsolator::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: response (gap effect)
  //------------------------------------------
  double KIsolator::gapResponse(double eta){
    double aEta = std::abs(eta);
    double p1 = 1.1095;
    double p2 = -7.494*std::pow((aEta - p1),2);
    double p3 = -25.275*std::pow((aEta - p1),2);
    double p4 = 14.52;
    return ((aEta>1.3)&&(aEta<1.9))? 
      1.0 - std::exp(p2 - p4*std::exp(p3)) : 1.0;
  }
  
  
  //-------------------------------------------
  // private: isIsolated
  //------------------------------------------
  
  // Checks whether a candidate electron is isolated from clusters
  // and cells
  
  bool KIsolator::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() 
                                      );
      
      // My own test: would be removed
      std::vector<ICluster*>::iterator it;
      for (it=clusters.begin(); it!=clusters.end(); ++it) {
        double dR = m_kinehelp.deltaR(*particle,*it);
        if (dR>m_rClusterMatch) continue;
        log << "\n\t -->  cluster: " << *it ;
        log << "\n\t -->   DeltaR: " << dR;
        log << "\n\t -->  EtC-EtP: " 
            << (*it)->eT() - 
          (*particle)->eT()*gapResponse((*particle)->eta());
        log << "\n\t ---------------------------------------- " <<endreq;
      }
      
      
      //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


      // dummy ReconstructedParticle to heve simple access to
      // unsmeared 4-vector to do isolation
      ReconstructedParticle* rp = 
        new ReconstructedParticle(
                                  (*particle)->pdg_id(),
                                  (*particle)->momentum(),
                                  (*particle)->truth()
                                  );


      // Use kinematic helper to do all work
      double eCoreSum = m_kinehelp.sumETInCone(
                                               storedCells.begin(), 
                                               storedCells.end(), 
                                               rp,
                                               m_rCellIsolation
                                               );
      

      // delete the dummy ReconstructedParticle
      delete rp;


      // 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() );


      // consider the effect of transition region (from FastShower)
      log << "\n\t -->Applying gap correction factor: "
          << gapResponse((*particle)->eta()) << endreq ; 

      double eCoreIsolation = eCoreSum;
      if(itDeposits) eCoreIsolation -=  
                       ((*particle)->eT())*gapResponse((*particle)->eta());
      
      
      // Apply isolation criteria
      
      if( eCoreIsolation  > m_eCellIsolation ) 
        log  
          << "\n\t -->Excess energy in core (R<0.2) = " << eCoreIsolation 
          << "   => NOT-ISOLATED " << endreq ;    
      else
        log   
          << "\n\t -->Excess energy in core (R<0.2) = " << eCoreIsolation
          << "   => ISOLATED " << endreq ;
                  
      if( eCoreIsolation  > 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
      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

---