JetMaker.cxx

Go to the documentation of this file.

// ================================================
// JetMaker class Implementation
// ================================================
//
// THIS TEXT TO BE REPLACED BY ATLAS STANDARD FORMAT
//
// Namespace Atlfast::
//
#include "AtlfastAlgs/JetMaker.h"
#include "AtlfastAlgs/ISmearer.h"
#include "AtlfastAlgs/JetSmearer.h"
#include "AtlfastAlgs/MissingMomentum.h"
#include "AtlfastAlgs/GlobalEventData.h"
#include "AtlfastAlgs/ContainerAssocsDispatcher.h"

#include "AtlfastEvent/Jet.h"
#include "AtlfastEvent/ICell.h"
#include "AtlfastEvent/CollectionDefs.h"
#include "AtlfastEvent/ReconstructedParticle.h"
#include "AtlfastEvent/ParticleCodes.h"
#include "AtlfastEvent/MsgStreamDefs.h"
#include "AtlfastEvent/ICluster.h"
#include "AtlfastEvent/MCparticleCollection.h"
#include "AtlfastEvent/SimpleKinematic.h"
#include "AtlfastEvent/Phi.h"

#include "AtlfastUtils/TesIO.h"
#include "AtlfastUtils/HeaderPrinter.h"
#include "AtlfastUtils/FunctionObjects.h"
#include "AtlfastUtils/IsAssociated.h"

#include <cmath> 
#include <algorithm>
#include <numeric>
#include <assert.h> // Gaudi includes
#include "GaudiKernel/DataSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"



// CLHEP,HepMC
#include "GeneratorObjects/McEventCollection.h"
#include "HepMC/GenEvent.h"
#include "HepMC/GenVertex.h"

#include "CLHEP/Units/SystemOfUnits.h"

namespace Atlfast{

  // STL helper function object    
    class SmearCell{
    public:
      SmearCell(ISmearer* is, double *sumET): 
        m_smearer(is), m_doSmear((is== 0)? false:true),m_sumET(sumET){}
      HepLorentzVector operator()(const HepLorentzVector& sumSoFar, 
                                  const ICell* c){
        
        HepLorentzVector sum(0., 0., 0., 0.);
        HepLorentzVector smeared;
        if(m_doSmear){
          smeared = m_smearer->smear( c->momentum() );
          sum = sumSoFar+smeared;
          (*m_sumET) += smeared.perp();
        }else{
          smeared = c->momentum(); // ie.. unsmeared
          sum = sumSoFar+smeared;
          (*m_sumET) += smeared.perp();
        }

        return sum; 
      }
    private:
      ISmearer* m_smearer;
      bool      m_doSmear;
      double   *m_sumET;
    };
      
    

//--------------------------------
// Constructors and destructors
//--------------------------------
  using std::abs;
  using std::sort;
  using std::partition;
  
  JetMaker::JetMaker 
  ( const std::string& name, ISvcLocator* pSvcLocator ) 
    : Algorithm( name, pSvcLocator ),
      m_bSelector(0),
      m_cSelector(0),
      m_tauSelector(0),
      m_smearer(0),
      m_cellSmearer(0),
      m_tesIO(0)
  {
    
    // Setting the parameter defaults.
    m_minPT            = 10*GeV;     // Min jet trans momentum
    m_maxEta           = 5.0;          //Max jet eta 
    m_doSmearing       = true;       //smearing is done
    m_rconeb           = 0.400;
    m_rconef           = 0.400;
    
    m_bPtMin           = 5.0*GeV;    //min pt for bquark to label jet
    m_bMaxDeltaR       = 0.3;        //max dR for bquark to label jet
    
    m_cPtMin           = 5.0*GeV;    //min pt for cquark to label jet
    m_cMaxDeltaR       = 0.3;
    
    m_tauPtMin         = 10.0*GeV;
    m_tauMaxDeltaR     = 0.3;
    
    m_etaTagMax        = 2.5;

    m_tauJetPtRatio    = 0;
    
    m_adjustMissETForIsolation = true;

    // Default paths for entities in the TES
    m_inputLocation             = "/Event/AtlfastClusters";
    m_outputLocation            = "/Event/AtlfastJets";
    m_muonLocation              = "/Event/AtlfastNonIsolatedMuons";
    m_unusedCellLocation        = "/Event/AtlfastUnusedCells";
    m_missingMomentumLocation   = "/Event/AtlfastMissingMomentum";
    m_isolatedElectronLocation  = "/Event/AtlfastIsolatedElectrons";
    m_isolatedPhotonLocation    = "/Event/AtlfastIsolatedPhotons";
    
    // This is how you declare the paramemters to Gaudi so that
    // they can be over-written via the job options file
    declareProperty( "MinimumPT",                m_minPT ) ;
    declareProperty( "MaximumEta",               m_maxEta ) ;
    declareProperty( "DoSmearing",               m_doSmearing ); 
    declareProperty( "RconeB",                   m_rconeb );
    declareProperty( "RconeF",                   m_rconef );
    declareProperty( "bPtMin",                   m_bPtMin );
    declareProperty( "bMaxDeltaR",               m_bMaxDeltaR );
    declareProperty( "cPtMin",                   m_cPtMin );
    declareProperty( "cMaxDeltaR",               m_cMaxDeltaR );
    declareProperty( "tauPtMin",                 m_tauPtMin );
    declareProperty( "tauMaxDeltaR",             m_tauMaxDeltaR );
    declareProperty( "etaTagMax",                m_etaTagMax );
    declareProperty( "tauJetPtRatio",            m_tauJetPtRatio);
    declareProperty( "adjustMissETForIsolation", m_adjustMissETForIsolation);

    declareProperty( "InputLocation",            m_inputLocation ) ;
    declareProperty( "OutputLocation",           m_outputLocation ) ;
    declareProperty( "UnusedCellLocation",       m_unusedCellLocation ) ;
    declareProperty( "MuonLocation",             m_muonLocation );
    declareProperty( "MissingMomentumLocation",  m_missingMomentumLocation );
    declareProperty( "IsolatedElectronLocation", m_isolatedElectronLocation );
    declareProperty( "IsolatedPhotonLocation",   m_isolatedPhotonLocation );

  }
  
  JetMaker::~JetMaker() {
    
    MsgStream log( messageService(), name() ) ;
    log << MSG::INFO << "Destructor Called" << endreq;
    
    if (m_smearer) {
      log << MSG::INFO << "Deleting jet Smearer" << endreq;
      delete m_smearer;
    }
    if (m_cellSmearer) {
      log << MSG::INFO << "Deleting Cell Smearer" << endreq;
      delete m_cellSmearer;
    }
    if (m_tesIO) {
      delete m_tesIO;
    }
  } 
  
  //---------------------------------
  // initialise() 
  //---------------------------------
  
  StatusCode JetMaker::initialize(){
    MsgStream log( messageService(), name() ) ;
    
    log << MSG::DEBUG << "instantiating a JetSmearer" << endreq;
    
    
    //get the Global Event Data using singleton pattern
    GlobalEventData* ged = GlobalEventData::Instance();
    int lumi                  = ged->lumi();
    int randSeed              = ged->randSeed() ;
    m_barrelForwardEta        = ged->barrelForwardEta();
    m_mcLocation              = ged->mcLocation();

    // Pass this into GlobalEventData for EventHeaderMaker to read later
    ged->set_adjustMissEtForIsolation(m_adjustMissETForIsolation);

    //    m_tesIO = new TesIO(eventDataService());
    m_tesIO = new TesIO(m_mcLocation, ged->justHardScatter());
    
    if(m_doSmearing){
      m_smearer                 = new JetSmearer(randSeed, 
                                                 lumi,
                                                 m_rconeb, 
                                                 m_rconef, 
                                                 m_barrelForwardEta);    
      m_cellSmearer             = new JetSmearer(randSeed, 
                                                 lumi,
                                                 0., 
                                                 0., 
                                                 m_barrelForwardEta);    
    } else {
      m_smearer = NULL;
      m_cellSmearer = NULL;
    }

    const int BQUARK(ParticleCodes::BQUARK);
    const int CQUARK(ParticleCodes::CQUARK);

    //.........................................................
    // Construct objects to select b, c and tau quarks with kine conditions
    
    HepMC_helper::IMCselector* selector;
    
    std::vector<HepMC_helper::IMCselector*> bSelectors;
    selector = new HepMC_helper::IsFromHardScatter();
    bSelectors.push_back( selector );
    selector = new HepMC_helper::SelectJetTag(BQUARK, m_bPtMin, m_etaTagMax);
    bSelectors.push_back( selector );
    m_bSelector = new HepMC_helper::NCutter(bSelectors);

    std::vector<HepMC_helper::IMCselector*> cSelectors;
    selector = new HepMC_helper::IsFromHardScatter();
    cSelectors.push_back( selector );
    selector = new HepMC_helper::SelectJetTag(CQUARK, m_cPtMin, m_etaTagMax);
    cSelectors.push_back( selector );
    m_cSelector = new HepMC_helper::NCutter(cSelectors);

    std::vector<HepMC_helper::IMCselector*> tauSelectors;
    selector = new HepMC_helper::IsFromHardScatter();
    tauSelectors.push_back( selector );
    selector = new HepMC_helper::SelectTauTag(m_tauPtMin, m_etaTagMax);
    tauSelectors.push_back( selector );
    m_tauSelector = new HepMC_helper::NCutter(tauSelectors);

    //.........................................................
    // Tidy up selectors
    std::vector<HepMC_helper::IMCselector*>::iterator iter;
    for(iter=bSelectors.begin();   iter!=bSelectors.end();   delete *iter, ++iter);
    for(iter=cSelectors.begin();   iter!=cSelectors.end();   delete *iter, ++iter);
    for(iter=tauSelectors.begin(); iter!=tauSelectors.end(); delete *iter, ++iter);
    
    

    HeaderPrinter hp("Atlfast JetMaker:", log);
    
    hp.add("TES Locations:              ");
    hp.add(" Muons from                 ", m_muonLocation); 
    hp.add(" Clusters from              ", m_inputLocation); 
    hp.add(" MC Input Location          ", m_mcLocation);
    hp.add(" Jets to                    ", m_outputLocation);
    hp.add(" unused cell,cluster sum to ", m_missingMomentumLocation);
    hp.add("Smearing on                 ", m_doSmearing);
    hp.add("Cluster min Pt              ", m_minPT);
    hp.add("Cluster max Eta             ", m_maxEta);           
    hp.add("Luminosity                  ", lumi);
    hp.add("Initial random seed         ", randSeed);      
    hp.add("End cap Cone size           ", m_rconef);
    hp.add("Barrel  Cone size           ", m_rconeb);
    hp.add("Eta of start of end cap     ", m_barrelForwardEta);
    hp.add("Parameters for labels       ", "       ");
    hp.add(" Max eta                    ", m_etaTagMax);
    hp.add(" b-quark min pt             ", m_bPtMin);
    hp.add(" b-quark max DeltaR         ", m_bMaxDeltaR);
    hp.add(" c-quark min pt             ", m_cPtMin);
    hp.add(" c-quark max DeltaR         ", m_cMaxDeltaR);
    hp.add(" tau lep min pt             ", m_tauPtMin);
    hp.add(" tau lep max DeltaR         ", m_tauMaxDeltaR);
    hp.add(" tau/Jet Pt Ratio           ", m_tauJetPtRatio);
    hp.print();
    


    return StatusCode::SUCCESS ;
  }
  
  //---------------------------------
  // finalise() 
  //---------------------------------
  
  StatusCode JetMaker::finalize(){

    MsgStream log( messageService(), name() ) ;
    log << MSG::INFO << "finalizing" << endreq;    

    if (m_bSelector)    delete m_bSelector;
    if (m_cSelector)    delete m_cSelector;
    if (m_tauSelector)  delete m_tauSelector;

    return StatusCode::SUCCESS ;
  }
  
  
  //----------------------------------------------
  // execute() method called once per event
  //----------------------------------------------
  //
  //  This algorithm creates smeared Jets passing cuts. 
  //  It scans the list of Clusters from the Monte Carlo truth. 
  //  If a Cluster of the right flavour is found it creates 
  //  a candidate Jet object by smearing the four vector
  //  of the Cluster and adding non-isolated muons. 
  // 
  //  It then applies kinematic criteria on the properties of the Jet
  //  Those which pass the cuts are kept.
  //  Finally all successful Jets are added to the event store.
  //
  //    //
  
  StatusCode JetMaker::execute( ){
    
    MsgStream log( messageService(), name() ) ;    

    std::string mess;
    // make a message logging stream
    
    //.............................
    // Make some locals stores which be used to keep pointers to
    // all of the entities from the event store which are needed by
    // this algorithm. These are all local and are typedefed in this class
    
    localClusterCollection      my_Clusters ;
    MuonCollection              my_Muons ;
    //................................
    
    //    if( ! m_tesIO->copy( my_Clusters, m_inputLocation ) ) {
    if( ! m_tesIO->copy<IClusterCollection> ( my_Clusters, m_inputLocation ) ) {
      log << MSG::INFO << "No Clusters in TES " << endreq;
    } else {
      log << MSG::DEBUG << my_Clusters.size()<<" Clusters from TES " << endreq;
    }
    // get non-isolated muons for jet construction
    
    if( ! m_tesIO->copy<ReconstructedParticleCollection>
        ( my_Muons, m_muonLocation ) ) {
      log << MSG::DEBUG << "No muons in TES " << endreq;
    }  

    // ......................
    // Make a container object which will store pointers to all isolated 
    // Jetss which are successfully created.  Since it is going to go 
    // into the event store, it has to be a special Gaudi type of container. 
    // This is typedefined in the entity include file
    // As it is going into the TES it needs to be created on the heap with new.
    
    JetCollection* myJets = new JetCollection ;
    
    //.........................................................
    // From each Cluster create a reconstructed Jet. 
    // If all requirements are satisfied add to the collection
    
    Jet* candidate ;
    localClusterIterator src ;
    m_sumET = 0.;
    HepLorentzVector missingMomentum(0.,0.,0.,0.);
    
    MCparticleCollection  mcParticles_b ;
    MCparticleCollection  mcParticles_c ;
    MCparticleCollection  mcParticles_tau ;

    //.........................................................
    // Get quarks used in tagging from TES using tesio

    TesIoStat stat = m_tesIO->getMC( mcParticles_b, m_bSelector ) ;
    if(!stat){
      log << MSG::ERROR <<"Problem getting b-quarks" << endreq; 
      return stat;
    }
    stat = m_tesIO->getMC( mcParticles_c, m_cSelector ) ;
    if(!stat){
      log << MSG::ERROR <<"Problem getting c-quarks" << endreq; 
      return stat;
    }
    stat = m_tesIO->getMC( mcParticles_tau, m_tauSelector ) ;
    if(!stat){
      log << MSG::ERROR <<"Problem getting taus" << endreq; 
      return stat;
    }
 
    for( src = my_Clusters.begin() ; src != my_Clusters.end() ; ++src ) { 
      //
      //      ReconstructedParticle rp;
      //      std::vector<IAOO*> v;
      //      v.push_back(*src);
      //      AssocTypeRecoverer atr(v.begin(), v.end());
      //      if ((atr.typeVector(rp))->size() == 0){
      if ( !IsAssociated<ReconstructedParticle>( *src ) ){
        // Make a new jet
        candidate = this->create( *src );
        
        //Make a temp HepLorentzVector from candidates 4-mom
        //before muon is added in
        HepLorentzVector temp4Vec = candidate->momentum();
        
        // add isolated muons
        this->addMuons(candidate, my_Muons);
        
        
        // add to missing momentum vector if not used in a jet.
        // remember, these are smeared momenta if smearing is enabled
        if( this->isAcceptable( candidate ) ) {    
          
          // if the jet is acceptable label the jet
          this->label(candidate,mcParticles_b,mcParticles_c,mcParticles_tau);
          
          log << MSG::DEBUG << "Pushing back a Jet" << endreq;
          myJets->push_back( candidate ) ; 
          log<<MSG::DEBUG<<"jet "<<*candidate<<endreq;
        }else {
          log << MSG::DEBUG << "Jet fails cuts" << endreq;
          //add in temp4Vec to avoid muon double counting
          missingMomentum += temp4Vec;
          m_sumET += temp4Vec.perp();
          delete candidate;
        }
      }else{
        log << MSG::DEBUG << "Cluster is associated, so ignored" << endreq;
      }
    }
    
    log<<MSG::DEBUG<<"Missing momentum: clusters "<<missingMomentum<<endreq;
    missingMomentum += addCells(log);
    log << MSG::DEBUG << "sumET after adding cells: " << m_sumET << endreq;
    log<<MSG::DEBUG<<"Missing momentum: cells+clstrs "
       <<missingMomentum<<endreq;

    //......................................
    // Register any Jets which have been successfilly created in the 
    // transient event store. Allow modification further along the line
    // for tagging routines.
    log<<MSG::DEBUG<<"Storing "<< myJets->size() <<" Jets"<<endreq;
    stat = m_tesIO->store( myJets, m_outputLocation, true) ;
    if(!stat){
      log<<MSG::ERROR<<"Could not store jets"<<endreq;
      return stat;
    }


    MissingMomentum* mm = new MissingMomentum(missingMomentum,m_sumET);
    stat = m_tesIO->store(mm,m_missingMomentumLocation ) ;
    std::string message = stat? "Stored missing mom":"Failed to store missing mom";
    log<<MSG::DEBUG<<message<<endreq;
  

    return StatusCode::SUCCESS;
  }
  
  
  
  //----------------------------------
  // create()
  //----------------------------------
  
  // Takes a single Cluster and creates a Jet
  // according to the desired criteria
  //
  // Note that we must NEW these jets so they can go to the TES
  // and if we decide that we do not want them, we must DELETE them
  // Once they are put to the TES, they are no longer our responsibility to delete
  
  Jet* JetMaker::create ( ICluster* src )
  {
    Jet* jet;
    MsgStream log( messageService(), name() ) ;
    HepLorentzVector vec = src->momentum();
    
    if (m_doSmearing) 
      {
        // Ask our Smearer make a smeared 4-vector
        vec = m_smearer->smear(vec);      
      }
    
    jet = new Jet(vec, *src);
    log << MSG::DEBUG << "Unsmeared Cluster: " << *src << endreq;
    log << MSG::DEBUG << "Smeared Jet      : " << jet << endreq ;
    
    //IAOO* js = jet;
    //assert(!js->unAssociated());
    // return Jet
    return jet;
    
    
  }
  //-------------------------------------------
  // addmuons
  //-------------------------------------------
  
  // adds non-isolated muons to the jets
  
  void JetMaker::addMuons(Jet* jet, MuonCollection& muons )
  {
    
    MsgStream log( messageService(), name() ) ;
    
    if (muons.size() == 0) return;
    
    MuonIterator first = muons.begin() ;
    MuonIterator last  = muons.end() ;
    
    // Partition muons to get those that satisfy the dR condition in either eta region
    MuonIterator partitionLast = 
      partition(first,last,PartitionCondition::BelowThresholdDeltaR(jet,m_rconeb,m_rconef,m_barrelForwardEta));
    
    if (partitionLast == first){
      log << MSG::DEBUG << "No muons in R cone" << endreq;
      return;
    }

    log << MSG::DEBUG << "Adding " << partitionLast-first << " muons to Jet" << endreq;

    // Add muon momenta to jet
    for (MuonIterator muItr = first; muItr != partitionLast;++muItr) jet->addMuon(*muItr);
    
    // Remove muons from list
    muons.erase(first,partitionLast);
    
  }
  
  
  //-------------------------------------------
  // isAcceptable
  //------------------------------------------
  
  // Decides whether a candidate is acceptable to this Maker, i.e. whether
  // it wants to proceed and add it to its output product collection
  
  bool JetMaker::isAcceptable ( Jet* candidate)
  {
    
    MsgStream log( messageService(), name() ) ;
    
    // Apply kinematic criteria to the candidate jet
    
    if( candidate->pT() < m_minPT ) {        
      log << MSG::DEBUG << "Candidate failed pt cut: pt " << candidate->pT() 
          << " min was " << m_minPT << endreq;
      return false ;
    }
    
    if( abs(candidate->eta()) > m_maxEta ) {
      log << MSG::DEBUG << "Candidate failed max eta cut: eta " 
          << candidate->eta() 
          << " max was " << m_maxEta << endreq;
      return false ;
    }
    
    return true ;
  }
  
  //-------------------------------------------
  // label
  //------------------------------------------
  //
  
  void JetMaker::label(Jet* jet, 
                       MCparticleCollection& mcParticles_b, 
                       MCparticleCollection& mcParticles_c, 
                       MCparticleCollection& mcParticles_tau)const {
    
    MsgStream log( messageService(), name() ) ;
    
    double dR = 0.;

    // The getJetLabel and getTauLabel functions loop over the input particle collections
    // and return the PDGID of a particle if within maxdeltaR of the jet (0 otherwise). 
    // dR is also set to the lowest found particle-jet delta R value

    int bLabel = getJetLabel(mcParticles_b, jet, m_bMaxDeltaR, dR);
    jet->setdRbquark(dR);
    int cLabel = getJetLabel(mcParticles_c, jet, m_cMaxDeltaR, dR);
    jet->setdRcquark(dR);
    int tauLabel = getTauLabel( mcParticles_tau, jet, m_tauMaxDeltaR, dR);
    jet->setdRhadtau(dR);

    // Set jet label (order dependent)
    if(bLabel){
      jet->setPdg_id(bLabel);
      return;
    }
    if(cLabel){
      jet->setPdg_id(cLabel);
      return;
    }
    if(tauLabel){
      jet->setPdg_id(tauLabel);
      return;
    }

    jet->setPdg_id(0);
    return; 

  }
  //-------------------------------------------
  // getJetLabel
  //------------------------------------------
  int JetMaker::getJetLabel(MCparticleCollection&  mcParticles, 
                            Jet* jet,
                            double maxDeltaR,
                            double &deltaR) const {
    
    MsgStream log( messageService(), name() ) ;

    int pdgId    = 0;
    double drMin = 999.;
    
    //Construct an object to select IKinematics within DeltaR of the current jet
    PartitionCondition::BelowThresholdDeltaR labelSelector(jet, maxDeltaR); 
    MCparticleCollectionCIter src ;
    for( src = mcParticles.begin(); src != mcParticles.end(); ++src ){
      SimpleKinematic sk(*src);
      if(labelSelector(sk)){
        pdgId = (*src)->pdg_id();
      }      
      Phi dPhi( sk.phi() - jet->phi() );
      double dR = sqrt((dPhi*dPhi) + 
                       (sk.eta()-jet->eta())*(sk.eta()-jet->eta()));
      if (dR < drMin) drMin = dR;
    }
    deltaR = drMin;
    return pdgId;
  }
  //-------------------------------------------
  // getTauLabel
  //------------------------------------------
  int JetMaker::getTauLabel(MCparticleCollection&  mcParticles,
                            Jet* jet,
                            double maxDeltaR,
                            double &deltaR) const {
    
    MsgStream log( messageService(), name() ) ;
    
    int pdgId    = 0;
    double drMin = 999.;

    //Construct an object to select IKinematics within DeltaR of the current jet
    PartitionCondition::BelowThresholdDeltaR labelSelector(jet, maxDeltaR);
    MCparticleCollectionCIter src ;
    for( src = mcParticles.begin() ; src != mcParticles.end() ; ++src )
      {
        HepLorentzVector tauHlv = (*src)->momentum();
        HepLorentzVector neuHlv;

        // Find Tau children
        HepMC::GenVertex*  endVertex = (*src)->end_vertex();
        if(!endVertex) continue;

        HepMC::GenVertex::particle_iterator firstChild = 
          endVertex->particles_begin(HepMC::children);
        
        HepMC::GenVertex::particle_iterator endChild = 
          endVertex->particles_end(HepMC::children);
        
        HepMC::GenVertex::particle_iterator thisChild = firstChild; 
        
        // Find Pt of Taus tau-neutrino
        for(; thisChild!=endChild; ++thisChild){
          if(  abs((*thisChild)->pdg_id()) == 16){
            neuHlv = (*thisChild)->momentum();
          }
        }
        double sigma;
        double sqrtene = sqrt(jet->momentum().e()/GeV);
        if ( fabs(jet->eta()) < m_barrelForwardEta )
          sigma = 0.5/sqrtene;
        else
          sigma = 1.0/sqrtene;

        double cutQuantity = (m_tauJetPtRatio) ? m_tauJetPtRatio : 1.0 - 2*sigma;

        log << MSG::DEBUG << "m_tauJetPtRatio = " << m_tauJetPtRatio
            << ", cutQuantity = " << cutQuantity << endreq;

        HepMC::GenParticle trueHadSystem;
        trueHadSystem.set_momentum(tauHlv-neuHlv); 
        SimpleKinematic sk(trueHadSystem);
        if (labelSelector(sk) &&
            (tauHlv-neuHlv).perp()/jet->pT() >= cutQuantity){
          Phi dPhi( sk.phi() - jet->phi() );
          deltaR = sqrt((dPhi*dPhi) +
                        (sk.eta()-jet->eta())*(sk.eta()-jet->eta()));
          pdgId = (*src)->pdg_id();
        }
        Phi dPhi( sk.phi() - jet->phi() );
        double dR = sqrt((dPhi*dPhi) + 
                         (sk.eta()-jet->eta())*(sk.eta()-jet->eta()));
        if (dR < drMin) drMin = dR;
      }
    
    deltaR = drMin;
    return pdgId;
  }   
  
  HepLorentzVector JetMaker::addCells(MsgStream& log){
    
    std::vector<const ITwoCptCell*> cellsToSmear(0);
    std::vector<const ITwoCptCell*> newCellsVector(0);
    
    getUnusedCells(cellsToSmear,log);
    
    if (m_adjustMissETForIsolation){
      getIsolationCorrectionCells(cellsToSmear,newCellsVector,m_isolatedElectronLocation,log);
      getIsolationCorrectionCells(cellsToSmear,newCellsVector,m_isolatedPhotonLocation,log);
    }
    
    log << MSG::DEBUG << "Smearing a total of " << cellsToSmear.size() << " cells" << endreq;
    log << MSG::DEBUG << "sumET before adding cells: " << m_sumET << endreq;
    
    // accumulate used to compute vectorial sum of Cell HepLorentzVectors
    // m_sumET also passed into SmearCell to compute sum of smeared cell ETs
    HepLorentzVector temp(0., 0., 0., 0.);
    HepLorentzVector sum = std::accumulate(cellsToSmear.begin(),
                                           cellsToSmear.end(),
                                           temp,
                                           SmearCell(m_cellSmearer,&m_sumET));
    //if (eHitCellPtr) delete eHitCellPtr;

    std::vector<const ITwoCptCell*>::iterator itccItr   = newCellsVector.begin();
    std::vector<const ITwoCptCell*>::iterator itccItrE  = newCellsVector.end();

    for (;itccItr!=itccItrE;++itccItr) delete *itccItr;
    newCellsVector.clear();

    return sum;
  }
  
  void JetMaker::getUnusedCells(std::vector<const ITwoCptCell*> &cellsToSmear, MsgStream& log){
    
    std::vector<const TwoCptCell*> unusedCells(0);
    if( !m_tesIO->copy<std::vector<const TwoCptCell*> >(unusedCells,
                                                        m_unusedCellLocation)){
      log << MSG::INFO << "No unused cells  in TES " << endreq;
      return;
    }
    
    log << MSG::DEBUG << unusedCells.size() << " unused cells in TES" << endreq;
    
    std::copy(unusedCells.begin(),unusedCells.end(),std::back_inserter(cellsToSmear));
    
    log << MSG::DEBUG  << "Added " << unusedCells.size() << " unused cells to cellToSmear" << endreq;
    
    return;
    
  }
  
  void JetMaker::getIsolationCorrectionCells(std::vector<const ITwoCptCell*> &cellsToSmear, 
                                             std::vector<const ITwoCptCell*> &newcellsVector,
                                             std::string tesAddress,
                                             MsgStream& log){
    
    
    std::vector<const ITwoCptCell*> isolationCorrectionCells(0);
    
    const ReconstructedParticleCollection* dh(0);
    if( !m_tesIO->getDH(dh, tesAddress) ) throw "Error in numberInList";
    
    ReconstructedParticleCollection::const_iterator iter = dh->begin();
    
    for (; iter != dh->end(); ++iter){
      
      const HepMC::GenParticle* gp = (*iter)->truth();

      // Get associations of ReconstructedParticle associations 
      const IAOO* iaooRP = (*iter);
      TypeVisitor tv =
        ContainerAssocsDispatcher( iaooRP->begin(), iaooRP->end(), TypeVisitor() );
      
      // Select the TwoCptCells
      std::vector<const TwoCptCell*> tccv = tv.typeVector(TwoCptCell());
      log << MSG::DEBUG << "TypeVisitor found " << tccv.size() 
          << " TwoCptCells in this cluster" << endreq;
      
      std::vector<const TwoCptCell*>::iterator tccvi;

      for ( tccvi = tccv.begin(); tccvi != tccv.end();){
        const ICell* ic = (*tccvi);

        // See if this TwoCptCell contains the GenParticle that made the
        // original ReconstructedParticle
        std::vector<const GenParticle*> gpv = ic->particles();
        std::vector<const GenParticle*>::iterator gpi;
        gpi = std::find(gpv.begin(),gpv.end(),gp);

        if ( gpi == gpv.end() ) ++tccvi;     // If not, go onto next TwoCptCell
        else{                                // It's in this one
          log << MSG::DEBUG << "GenParticle found in this cell, one of " 
              << gpv.size() << " GenParticles" << endreq;
          log << MSG::DEBUG << "GenParticle momentum = " << (*gpi)->momentum() << endreq;
          log << MSG::DEBUG << "Cell momentum = " << ic->momentum() << endreq;

          if ( gpv.size() > 1){
            // Make a new cell with pT reduced by that of the originating GenParticle
            ITwoCptCell* newcell_itcc = (*tccvi)->cloneITCC();
            newcell_itcc->resetCell();
            ICell* newcell = newcell_itcc;
            double newpT = ic->momentum().perp() - (*gpi)->momentum().perp();
            newpT = ( fabs(newpT) < 1e-10 ) ? 0 : newpT; // problem with precision
            newcell->newHit( newpT );
            
            // Add the new, adjusted hit cell
            isolationCorrectionCells.push_back(newcell_itcc);
            newcellsVector.push_back(newcell_itcc);
          }
          
          // Remove the pointer to the hit cell
          tccvi = tccv.erase(tccvi);
          
        }
      }
      // After adjustment, copy all cells into isolationCorrectionCells
      std::copy(tccv.begin(),tccv.end(),std::back_inserter(isolationCorrectionCells));
    }
    
    log << MSG::DEBUG << "Added " << isolationCorrectionCells.size() 
        << " isolation correction cells to cellsToSmear" << endreq;
   
    // Add isolationCorrectionCells cells to vector of cells for smearing
    std::copy(isolationCorrectionCells.begin(),
              isolationCorrectionCells.end(),
              std::back_inserter(cellsToSmear));
    
    return;
    
  }
  
} // end namespace bracket

---