// ================================================
 // Isolator class Implementation
 // ================================================
 //
 // THIS TEXT TO BE REPLACED BY ATLAS STANDARD FORMAT
 //
 // Namespace ATLFast
 //
 
 #include "AtlfastAlgs/Isolator.h"
 #include "AtlfastAlgs/JetSmearer.h"
 #include "AtlfastAlgs/GlobalEventData.h"
 
 #include "AtlfastEvent/Phi.h"
 #include "AtlfastEvent/ICluster.h"
 
 #include "AtlfastUtils/IsAssociated.h"
 #include "AtlfastUtils/FunctionObjects.h"
 #include "AtlfastUtils/HeaderPrinter.h"
 #include "TruthHelper/GenIMCselector.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;
 
 /**
    * @brief Associates ReconstructedParticles to Clusters and evaluates isolation.
    *
    * The associations are established by searching in a user-defined cone around
    * the RP position for unassociated Clusters.<br>
    * The isolation is calculated by searching in user-defined cones around the
    * RP position for Clusters and unused (unclustered) Cells. The corresponding
    * energies are then compared to user-defined limits.
    *
    * @image html Isolator.jpg "Isolator sequence"
    */
 
   //--------------------------------
   // Constructors and destructors
   //--------------------------------
   
   Isolator::Isolator 
   ( const std::string& name, ISvcLocator* pSvcLocator ) 
     : Algorithm( name, pSvcLocator ),
       m_cellSmearer( 0 ),
       m_tesIO( 0 ) {
     
     // Setting the parameter defaults.
     /** R-cone for associating a Cluster to a particle */
     m_rClusterMatch      = 0.150;
     /** R-cone in which to sum unassociated Cluster eT */
     m_rClusterIsolation  = 0.400;
     /** eT threshold to apply to summed Cluster eT */
     m_eClusterIsolation  = 0.0*GeV;
     /** R-cone in which to sum Cell eT */
     m_rCellIsolation     = 0.200;
     /** eT threshold to apply to excess summed Cluster eT */
     m_eCellIsolation     = 10.0*GeV;
     
     /** Upper and lower bounds of halo interval in dR */
     m_rLowerHalo         = 0.2;
     m_rHigherHalo        = 0.4;
     
      /** Whether or not to smear cells in isolation ET sum */
     m_smearCells         = false;
 
     m_fastShower         = false;
 
     // Default paths for entities in the TES
 
     /** Location in TES to obtain input particles to be isolated */
     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 ) ;
 
     declareProperty( "SmearCells",                m_smearCells ) ;
     declareProperty( "FastShower",                m_fastShower ) ;
   }
   
   // Destructor
   Isolator::~Isolator() {
     MsgStream log( messageService(), name() ) ;
     log << MSG::INFO << "destructor" << endreq;
     if (m_cellSmearer) {
       log << MSG::INFO << "Deleting Cell Smearer" << endreq;
       delete m_cellSmearer;
     }
     if (m_tesIO) {
       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();
     int randSeed             = ged->randSeed() ;
     int lumi                 = ged -> lumi();
     double barrelForwardEta  = ged -> barrelForwardEta();
 
     m_tesIO = new TesIO(m_mcLocation, ged->justHardScatter());
 
     if ( m_smearCells ){
       m_cellSmearer             = new JetSmearer(randSeed, 
                                                  lumi,
                                                  0., 
                                                  0., 
                                                  barrelForwardEta);    
     } else {
       m_cellSmearer = 0;
     }
     
     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.add("Smear Cells in isolation?  ", m_smearCells);
     hp.add("Running FastShower         ", m_fastShower);
     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 reconstructed particles (electron, photon or muon) 
   //  It scans the list of MC reconstructed particles. If a particle is found it creates 
   //  a candidate particle 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 reconstructed particle kept.
   //  Finally all successful reconstructed particles 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 particles.  
     
     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 Isolator::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 particle 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 ;
     
     // Measure summed cell ET within dR cones
     // To be used for analysis-level isolation cuts
     measureCones(log, 
                  particle, 
                  !storedCells.empty(),
                  storedCells.begin(), 
                  storedCells.end() );
     
     // 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 = 0 ;
     
 
 
     //....................................
     // 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() 
                                       );
      
       if(m_fastShower){  
          // 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() ) {
     if(!m_fastShower){  
          
       // 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 ;
       
      } else { // do the following for fast shower
       
       // 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 ));
    } // end of associated
      
     log << endreq ;
     
     return true ;
   } // end of isIsolated method
 
   void Isolator::measureCones(
                               MsgStream& log,
                               ReconstructedParticleCollection::const_iterator particle,
                               bool useCells,
                               std::vector<ITwoCptCell*>::iterator firstCell,
                               std::vector<ITwoCptCell*>::iterator lastCell
                               ){
     
     // Initialise etcone variables
     double etcone10 = 0, etcone20 = 0, etcone30 = 0, etcone40 = 0;
     
     //............................................
     // Cell isolation: 
     
     if ( useCells ){
       
       if ( m_smearCells && m_cellSmearer ){
         
         etcone10 += m_kinehelp.sumSmearedETInCone( firstCell, lastCell, *particle, 0.1, m_cellSmearer );
         etcone20 += m_kinehelp.sumSmearedETInCone( firstCell, lastCell, *particle, 0.2, m_cellSmearer );
         etcone30 += m_kinehelp.sumSmearedETInCone( firstCell, lastCell, *particle, 0.3, m_cellSmearer );
         etcone40 += m_kinehelp.sumSmearedETInCone( firstCell, lastCell, *particle, 0.4, m_cellSmearer );
         
 
       } else {
         etcone10 += m_kinehelp.sumETInCone( firstCell, lastCell, *particle, 0.1 );
         etcone20 += m_kinehelp.sumETInCone( firstCell, lastCell, *particle, 0.2 );
         etcone30 += m_kinehelp.sumETInCone( firstCell, lastCell, *particle, 0.3 );
         etcone40 += m_kinehelp.sumETInCone( firstCell, lastCell, *particle, 0.4 );
       }
       
       // 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){
         double particleET = (*particle)->eT();
         etcone10 -= particleET ;
         etcone20 -= particleET ;
         etcone30 -= particleET ;
         etcone40 -= particleET ;
       }
     
     }
     
     (*particle)->set_etcone10(etcone10);
     (*particle)->set_etcone20(etcone20);
     (*particle)->set_etcone30(etcone30);
     (*particle)->set_etcone40(etcone40);
     
     log << MSG::DEBUG << "etcone10 = " << (*particle)->etcone10()
         << ", etcone20 = " << (*particle)->etcone20()
         << ", etcone30 = " << (*particle)->etcone30()
         << ", etcone40 = " << (*particle)->etcone40()
         << endreq;
   }
   
 } // end of namespace bracket
 

Due to the LXR bug, the updates fail sometimes to remove references to deleted files. The Saturday's full rebuilds fix these problems
This page was automatically generated by the LXR engine.