TauMaker.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------
//
// file:     TauMaker.cxx
//
// Creates Tau candidates from Atlfast Tracks. 
//
//-----------------------------------------------------------------------
#include <fstream>
#include <cmath>

#include <GaudiKernel/MsgStream.h>

#include "CLHEP/Vector/LorentzVector.h"
#include "CLHEP/Random/JamesRandom.h"
#include "CLHEP/Random/RandGauss.h"
#include "CLHEP/Random/RandFlat.h"
#include "CLHEP/Units/SystemOfUnits.h"
#include "CLHEP/Random/RandomEngine.h"

#include "HepMC/GenVertex.h"

#include "AtlfastEvent/Phi.h"
#include "AtlfastEvent/Track.h"
#include "AtlfastEvent/CollectionDefs.h"

#include "AtlfastAlgs/TauMaker.h"
#include "AtlfastAlgs/GlobalEventData.h"

#include "PathResolver/PathResolver.h"

// -------------------------------------

namespace Atlfast
{
  //-------------------------------------------------------------------------
  // Constructor
  //-------------------------------------------------------------------------
  TauMaker::TauMaker( const std::string& name, ISvcLocator* pSvcLocator ) 
    : Algorithm( name, pSvcLocator ),
      m_tesIO(0)
  {
    
    //Default private parameters   
    m_TrackEffFile       = "atlfastDatafiles/trackEffInputFile.txt";
    m_TrackCollection    = "/Event/AtlfastTracks";
    m_outputLocation     = "/Event/AtlfastTaus";
    m_recoEtaCut         = 2.5; 
    m_pTLeadTrackCut     = 9000.0;
    m_pTOtherTrackCut    = 1000.0;
    m_maxNumTracks       = 3;
    m_detRCoreTrackCut   = 0.2;
    m_detRCoreCaloCut    = 0.2;
    m_detRRecoSeparation = 0.4; //Min separation of seed track from other reco tau
    m_detRTrueTauCut     = 0.2;
    m_recoETCut          = 10000.0;
    
    m_resEflowEne        = 0.3; //Energy Smearing: Need to optimise - guess 30% for now
    m_resNeuhEne         = 0.3; //Neutron ECal response - guess 30% for now...
    m_resTrack           = 0.5; //Keep for now, may remove later
    m_eleVetoEff         = 0.92; //Numbers adapted from TP Week June 2007
    m_muoVetoEff         = 0.96; //Numbers adapted from TP Week June 2007
    
    
    //Declare joboption properties here...
    declareProperty( "TrackCollection",    m_TrackCollection    );
    declareProperty( "OutputLocation",     m_outputLocation     );
    declareProperty( "RecoEtaCut",         m_recoEtaCut         );
    declareProperty( "pTLeadTrackCut",     m_pTLeadTrackCut     );
    declareProperty( "pTOtherTrackCut",    m_pTOtherTrackCut    );
    declareProperty( "maxNumTracks",       m_maxNumTracks       );
    declareProperty( "detRCoreTrackCut",   m_detRCoreTrackCut   );
    declareProperty( "detRCoreCaloCut",    m_detRCoreCaloCut    );
    declareProperty( "detRRecoSeparation", m_detRRecoSeparation );
    declareProperty( "detRTrueTauCut",     m_detRTrueTauCut     );
    declareProperty( "resEflowEne",        m_resEflowEne        );
    declareProperty( "recoETCut",          m_recoETCut          );
    declareProperty( "resNeuhEne",         m_resNeuhEne         );
    declareProperty( "resTrack",           m_resTrack           ); //Possibly obsolete
    declareProperty( "eleVetoEff",         m_eleVetoEff         );
    declareProperty( "muoVetoEff",         m_muoVetoEff         );
    declareProperty( "TrackEffFile",       m_TrackEffFile       );
  }
  
  TauMaker::~TauMaker() {
    
    if (m_tesIO) {
      delete m_tesIO;
    }
  }
  
  //----------------------------------------------------------------
  // Tool initialization
  //----------------------------------------------------------------
  StatusCode TauMaker::initialize()
  {
    MsgStream p_log( messageService(), name() ) ;

    p_log << MSG::VERBOSE << "RecoEtaCut         = ";
    p_log << m_recoEtaCut << endreq;
    
    
    GlobalEventData* ged = GlobalEventData::Instance();
    int randSeed = ged->randSeed() ;
    m_tesIO      = new TesIO(ged->mcLocation(), ged->justHardScatter());
    
    // Random Number Generators
    // Concerns about using the same random number 
    // engine for both m_randGauss and m_randFlat.
    // pass engine by reference rather than by pointer
    // to avoid deletion issues. (if passed by pointer
    // the engine is deleted when m_randGauss is deleted,
    // but then you can't delete m_randFlat without seg fault)
    // OR MAYBE USE GLOBAL RANDOM NUMBER ENGINE?
    
    m_randEngine   = new HepJamesRandom( randSeed );
    m_randGauss    = new RandGauss( *m_randEngine );
    m_randFlat     = new RandFlat( *m_randEngine );
    
    
    p_log << MSG::DEBUG << "m_TrackEffFile before PathResolver: " 
          << m_TrackEffFile <<endreq;
    m_TrackEffFile = PathResolver::find_file (m_TrackEffFile, "DATAPATH");
    p_log << MSG::DEBUG << "m_TrackEffFile after PathResolver: " 
          << m_TrackEffFile <<endreq;
    
    // Read in track reconstruction efficiency parametrisation
    std::ifstream trackEffFile;
    trackEffFile.open(m_TrackEffFile.c_str());
    
    if(trackEffFile){    
      p_log << MSG::INFO 
            << "Track Efficiency File " 
            << m_TrackEffFile
            << " open."
            << endreq;    
      double pTpoint  = 0.;
      double hadTrackEff = 0.;
      double muoTrackEff = 0.;
      double eleTrackEff = 0.;
      int ncolumn;        
      trackEffFile>>ncolumn;
      // Read corr factor bins from file
      for(int i=0; i<ncolumn; i++){
        trackEffFile>>pTpoint>>hadTrackEff>>muoTrackEff>>eleTrackEff;
        m_trackPTPoint.push_back(pTpoint);
        m_trackHadEff.push_back(hadTrackEff);
        m_trackMuoEff.push_back(muoTrackEff);
        m_trackEleEff.push_back(eleTrackEff);
      }
    }else{
      p_log << MSG::ERROR 
            << "Track Efficiency File file not found" 
            <<endreq;
      return StatusCode::FAILURE ;
    }
    trackEffFile.close();
    
    if ( m_trackPTPoint.size() != m_trackHadEff.size() ){
      p_log << MSG::FATAL;
      p_log << "Inconsistent Size of Track Efficiency Vectors";
      p_log << endreq;
    }
    
    return StatusCode::SUCCESS;
  }
  
  //---------------------------------------------------------------
  // Tool execution
  //---------------------------------------------------------------
  StatusCode TauMaker::execute()
  {

    MsgStream p_log( messageService(), name() ) ;
    
    // Container for new Tau objects
    TauContainer *container = new TauContainer;

    //---------------------------------------------------------
    // retrieve container with MC particles
    //---------------------------------------------------------
    std::vector<const HepMC::GenParticle*> p_particles;
    
    TesIoStat stat = m_tesIO->getMC( p_particles);
    
    if( !stat )
      {
        p_log << MSG::ERROR << "Unable to getMC";
        p_log << endreq;
        return StatusCode::RECOVERABLE;
        
      }
    
    //----------------------------------------------------------
    // retrieve pointer with Atlfast track collection
    //----------------------------------------------------------
    
    // get Atlfast track collection
    const TrackCollection* atlColl = 0;
    
    if(!m_tesIO->getDH(atlColl, m_TrackCollection)){
      p_log << MSG::DEBUG
            << "Could not retrieve at [" << m_TrackCollection << "] !!"
            << endreq ;
      return StatusCode::RECOVERABLE ;
    }
    
    if ( p_log.level() <= MSG::WARNING ) {
      p_log << MSG::DEBUG 
            << "Number of Atlfast tracks: [" 
            << atlColl->size() << "] ..."
            << endreq;
    }
    
    const TrackCollection::const_iterator itrEnd = atlColl->end();
    
    //-------------------------------------------------------------
    // Randomly tag tracks as having "good reconstruction"
    // pT dependent parametrisation of efficiency
    
    std::vector<unsigned int> TrackQualityVec;    
    // 1) First loop over tracks
    for ( TrackCollection::const_iterator itr_tag = atlColl->begin();
          itr_tag != itrEnd;
          ++itr_tag ) {
      
      const Track * atlTrackTag = *itr_tag;
      unsigned int qualitytag = 0;
      // 2) Draw random number    
      float rtag = randFlat()->fire();
      // 3) Get efficiency for good recon according to pT
      if( rtag < GetTrackEfficiency(atlTrackTag) ) { qualitytag = 1; }
      // 4) fill vector of tags identifying well reconstructed tracks
      TrackQualityVec.push_back( qualitytag );
    }

    // Check that the size of the TrackQualityVec vector matches
    // the size of the track collection.
    if( TrackQualityVec.size() != atlColl->size() ){
      p_log << MSG::WARNING << "TrackQualityVec Vector Not of Correct Size" <<endreq;     
    }

    //---------------------------------------------------------
    // Loop over Tracks in track collection
    //---------------------------------------------------------

    int seedTrackCounter = 0;
        
    for ( TrackCollection::const_iterator itr = atlColl->begin();
          itr != itrEnd;
          ++itr ) {
      int tr_IsTrueTau      = 0;

      double resphi         = 0; 

      // ET and position of tau candidate
      double tr_ET          = 0;
      double tr_Eta         = 0;
      double tr_Phi         = 0;

      // Eta Phi position of seed track
      double tr_SeedPhi     = 0;
      double tr_SeedEta     = 0;

      // Component From Tracks
      int tr_NTrackGood      = 0;
      int tr_NTrackPoor      = 0;
      int tr_SumQTrackGood   = 0;
      int tr_SumQTrackPoor   = 0;

      double tr_PTTrack      = 0;
      double tr_EtaTrack     = 0;
      double tr_PhiTrack     = 0;
          
      double tr_PTTrackMax   = 0;
      double tr_EtaTrackMax  = 0;
      double tr_PhiTrackMax  = 0;
            
      // Neutral Component of Energy
      int tr_NPi0           = 0;
      double tr_ETPi0       = 0;
      double tr_EtaPi0      = 0;
      double tr_PhiPi0      = 0;

      // True Leptonic Tracks...
      int tr_NMuoTrack = 0;
      int tr_NEleTrack = 0;
      // Lepton veto flags
      int tr_muoVetoFlag = 0;
      int tr_eleVetoFlag = 0;

      const Track * atlTrack = *itr;

      p_log << MSG::VERBOSE << "Search Seed Track "<<seedTrackCounter<<"\tquality = "
            <<TrackQualityVec.at(seedTrackCounter)<<"\tpdgid = "<<atlTrack->pdg_id()<<"\t"<<endreq;

      int seedTrackQuality = TrackQualityVec.at(seedTrackCounter);
      seedTrackCounter++;

      //-------------------------------------------
      // Select Primary Track from Track Collection
      //-------------------------------------------

      //-----------------------------------------------
      // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      // Temporary fix for Track bug in release 13.0.10
      // Ignore tracks with |PID|=24
      //-----------------------------------------------
      if( abs(atlTrack->pdg_id() )== 24 ) {
        continue; }

      //---------------------------------------------------------
      // remove tracks not tagged as having "good reconstruction"
      //---------------------------------------------------------
      if( seedTrackQuality == 0 ) {
        continue; }
            
      //-------------------------------------------
      // remove candidates with track momenta below the threshold
      //-------------------------------------------
      if(  atlTrack->pT() < m_pTLeadTrackCut  ) {
        continue; }

      //-------------------------------------------
      // remove candidates with pseudorapidity outside range
      //-------------------------------------------
      if(  std::fabs( atlTrack->eta() ) > m_recoEtaCut  ) {
        continue; }

      p_log << MSG::DEBUG << "---> Seed Track ( pT= "<<atlTrack->pT()/GeV
            <<" GeV, eta= "<<atlTrack->eta()
            <<", phi= "<<atlTrack->phi()
            <<", PID= "<<atlTrack->pdg_id()<<" )."<<endreq;

      //-------------------------------------------
      // look to see if it overlaps with candidaes
      // already in TauContainer
      //-------------------------------------------
      int separated = 1;

      TauContainer::const_iterator p_iterF = container->begin();
      for( ; p_iterF != container->end(); p_iterF++ )
        {
          //dR separation of seed track from other reconstructed taus
          Phi    p_detPhi( atlTrack->phi() - (*p_iterF)->phi() );
          double p_detEta = fabs( atlTrack->eta() - (*p_iterF)->eta() );
          double p_detCone = sqrt((p_detPhi * p_detPhi) + (p_detEta * p_detEta) );

          if( p_detCone <  m_detRRecoSeparation )separated = 0;
        }
          
      if( separated == 0 ) continue;


      //-------------------------------------------
      // Selected track accepted as a leading seed
      //-------------------------------------------
          
      tr_NTrackGood++;
      tr_SumQTrackGood = atlTrack->trajectory().signOfCharge(); 
      if( abs(atlTrack->pdg_id()) == 11 ) { tr_NEleTrack++; }
      if( abs(atlTrack->pdg_id()) == 13 ) { tr_NMuoTrack++; }

      tr_PTTrack       = atlTrack->pT();
      tr_EtaTrack      = atlTrack->eta() * atlTrack->pT();
          
      tr_ET  = atlTrack->eT();
      tr_Eta = atlTrack->eta() * atlTrack->pT();

      tr_SeedEta = atlTrack->eta();
      tr_SeedPhi = atlTrack->phi();
          
      tr_PTTrackMax  = atlTrack->pT();
      tr_EtaTrackMax = atlTrack->eta();
      tr_PhiTrackMax = atlTrack->phi();
          
      //-------------------------------------------
      // Start iterating over other tracks
      //-------------------------------------------
      int assocTrackCounter = 0;

      for ( TrackCollection::const_iterator itr2 = atlColl->begin();
            itr2 != itrEnd;
            ++itr2 ) {

        const Track * atlTrack2 = *itr2;
            
        //---------------------------------------------------------
        // does track have "good reconstruction"?
        //---------------------------------------------------------
        p_log << MSG::VERBOSE << "Associated Track "<<assocTrackCounter<<"\t \"quality\"="
              <<TrackQualityVec.at(assocTrackCounter)<<endreq;

        int assocTrackQuality = TrackQualityVec.at(assocTrackCounter);
        assocTrackCounter++;

        //-----------------------------------------------
        // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        // Temporary fix for Track bug in release 13.0.10
        // Ignore tracks with |PID|=24
        //-----------------------------------------------
        if( abs(atlTrack2->pdg_id() )== 24 ) {
          continue; }

        //---------------------------------------------------------
        // Check not double counting tracks
        //---------------------------------------------------------
        if( atlTrack2 == atlTrack ) {
          continue; }

        //---------------------------------------------------------
        // remove candidates with pseudorapidity outside range
        //---------------------------------------------------------
        if(  std::fabs( atlTrack2->eta() ) > m_recoEtaCut  ) {
          continue; }

        //dR separation of track from the seed track
        Phi    tr_detPhi( atlTrack->phi() - atlTrack2->phi() );
        double tr_detEta = fabs( atlTrack->eta() - atlTrack2->eta() );
        double tr_detCone = sqrt( (tr_detPhi * tr_detPhi) + (tr_detEta * tr_detEta) );

        p_log << MSG::VERBOSE << "Search Associated Track pT = "<<atlTrack2->pT()
              <<"\t dR = "<<tr_detCone<<endreq;

        if ( tr_detCone < m_detRCoreTrackCut ) {

          // Track Passes pT cut...
          if (atlTrack2->pT() > m_pTOtherTrackCut ){

            // Good Core Tracks:                
            if ( assocTrackQuality == 1 ){
              p_log << MSG::DEBUG << "-----> \"Good\" Track in Core Cone pT = "
                    <<atlTrack2->pT()/GeV << "\t detCone = "<<tr_detCone
                    <<"\tPID= "<< atlTrack2->pdg_id()<<endreq;

              tr_NTrackGood++;
              tr_SumQTrackGood += atlTrack2->trajectory().signOfCharge();
            }
            // Other Core tracks:
            if ( assocTrackQuality == 0 ){
              p_log << MSG::DEBUG << "-----> \"Poor\" Track in Core Cone pT = " 
                    << atlTrack2->pT()/GeV <<"  \tdetCone = "<<tr_detCone
                    <<" \tPID= "<< atlTrack2->pdg_id()<<endreq;
                  
              tr_NTrackPoor++;
              tr_SumQTrackPoor += atlTrack2->trajectory().signOfCharge();
            }

            tr_ET  += atlTrack2->eT();
            tr_Eta += atlTrack2->eta()* atlTrack2->pT();
            // resphi = difference in phi between othertrack and seed track 
            
            resphi  = Phi(atlTrack2->phi() - tr_SeedPhi);
            tr_Phi += Phi(resphi * atlTrack2->pT());

            if( abs(atlTrack2->pdg_id()) == 11 ) { tr_NEleTrack++; }
            if( abs(atlTrack2->pdg_id()) == 13 ) { tr_NMuoTrack++; }

            tr_PTTrack  += atlTrack2->pT();
            tr_EtaTrack += atlTrack2->eta()* atlTrack2->pT();
            resphi  = Phi(atlTrack2->phi() - tr_SeedPhi);
            tr_PhiTrack += resphi * atlTrack2->pT();            

            if( tr_PTTrackMax < atlTrack2->pT() )
              {
                tr_PTTrackMax  = atlTrack2->pT();
                tr_EtaTrackMax = atlTrack2->eta();
                tr_PhiTrackMax = atlTrack2->phi();
              }
          }
              
              
          // Low pT tracks in core region
          else if(atlTrack2->pT() <= m_pTOtherTrackCut){
            p_log << MSG::DEBUG << "-----> \"Low pT\" Track in Core Cone pT = " << atlTrack2->pT()/GeV 
                  <<" \tdetCone = "<<tr_detCone
                  <<" \tPID= "<< atlTrack2->pdg_id()<<endreq;
          }
        }
      } // end of itr2 loop

      p_log << MSG::DEBUG << "-----> NTrackGood = " <<tr_NTrackGood 
            << "\tPTTrack = " <<tr_PTTrack/GeV << endreq; 
      p_log << MSG::DEBUG << "-----> NTrackPoor = " <<tr_NTrackPoor << endreq;
      p_log << MSG::DEBUG << "-----> Of which NMuoTrack = "<<tr_NMuoTrack 
            << "; and NEleTrack = " <<tr_NEleTrack << endreq; 
              
      //--------------------------------------------------------------

      // m_maxNumTracks set to 3 by default. Can be changed in joboptions
      // to handle candidates with higher track multiplicities

      //if Ntracks is greater than maximum number of tracks cut, reject and continue
      if( tr_NTrackGood > m_maxNumTracks ) continue;
          
      p_log << MSG::DEBUG << "Good Tracks SumQ = " <<tr_SumQTrackGood 
            << "; Poor Tracks SumQ = " <<tr_SumQTrackPoor << endreq;

      // Sum Charge of tracks
      int tr_SumQ = tr_SumQTrackGood;
      if( tr_NTrackGood == 2 && tr_NTrackPoor==1 ) tr_SumQ = tr_SumQTrackGood + tr_SumQTrackPoor;

      int tr_NProng          = 0;
      //Decide whether its a 1-prong, 2-prong, 3-prong etc...
      if( tr_NTrackGood==1 )                     tr_NProng=1;
      if( tr_NTrackGood==2 && tr_NTrackPoor==0 ) tr_NProng=2;
      if( tr_NTrackGood==2 && tr_NTrackPoor==1 ) tr_NProng=3;
      if( tr_NTrackGood==2 && tr_NTrackPoor>=2 ) continue;
      if( tr_NTrackGood==3 )                     tr_NProng=3;
      //and for higher n track candidates...
      if( tr_NTrackGood >3 && tr_NTrackGood<=m_maxNumTracks) tr_NProng=tr_NTrackGood;

      // Reject Candidates with bad charge sum
      if( tr_NProng == 3 && abs(tr_SumQ) != 1 ) 
        {
          p_log << MSG::DEBUG << "Candidate Rejected NProng =" << tr_NProng
                << "  SumQ = " << tr_SumQ << endreq;
          continue;           
        }

      //--------------------------------------------------------------
      // Having found associated tracks, now look for neutral content
      // loop over truth container and find nearby pions
      //--------------------------------------------------------------

      std::vector<const HepMC::GenParticle*>::iterator p_MC_Iter;
          
      for( p_MC_Iter = p_particles.begin(); p_MC_Iter != p_particles.end();
           p_MC_Iter++ ) {
        if( (*p_MC_Iter)->barcode() >= 99999. ) {
          continue; }
            
        //-------------------------------------------------------
        // Check that particle is inside eta cut??
        //-------------------------------------------------------
        if(  std::fabs(  (*p_MC_Iter)->momentum().pseudoRapidity()) > m_recoEtaCut  )
          continue;
            
        //-------------------------------------------------------
        // Find separation of neutral pion from seed track
        //-------------------------------------------------------
        double p_detPhi = Phi( atlTrack->phi() - (*p_MC_Iter)->momentum().phi() );
        double p_detEta = fabs( atlTrack->eta() - (*p_MC_Iter)->momentum().pseudoRapidity() );
        double p_detCone = sqrt((p_detPhi * p_detPhi) + (p_detEta * p_detEta) );            

        //-------------------------------------------------------
        // neutral pions content
        //-------------------------------------------------------
        if( std::abs( (*p_MC_Iter)->pdg_id() ) == 111 ){

          p_log << MSG::VERBOSE << "Pi0 pT= " << (*p_MC_Iter)->momentum().perp()/GeV 
                << "\teta = " << (*p_MC_Iter)->momentum().eta() 
                <<"\tphi= "<<(*p_MC_Iter)->momentum().phi()<<"\tdR= "<<p_detCone<< endreq;
              
          if(  p_detCone <  m_detRCoreCaloCut ) {

            p_log << MSG::DEBUG << "---> Core Pi0 pT= " << (*p_MC_Iter)->momentum().perp()/GeV 
                  << "\teta = " << (*p_MC_Iter)->momentum().eta() <<"\tphi= "
                  <<(*p_MC_Iter)->momentum().phi()<< endreq;

            tr_ET  += (*p_MC_Iter)->momentum().perp();
            tr_Eta += (*p_MC_Iter)->momentum().pseudoRapidity() * (*p_MC_Iter)->momentum().perp();
            resphi  = Phi( (*p_MC_Iter)->momentum().phi() - tr_SeedPhi );
            tr_Phi += resphi * (*p_MC_Iter)->momentum().perp();

            tr_NPi0++;

            tr_ETPi0  += (*p_MC_Iter)->momentum().perp();
            tr_EtaPi0 += (*p_MC_Iter)->momentum().pseudoRapidity() * (*p_MC_Iter)->momentum().perp();
            resphi  = Phi( (*p_MC_Iter)->momentum().phi() - tr_SeedPhi );
            tr_PhiPi0 += resphi * (*p_MC_Iter)->momentum().perp();
                
          }
        } else if( std:: abs((*p_MC_Iter)->pdg_id()) == 2112 ){
          // Energy from neutrons in core cone...
          if(  p_detCone <  m_detRCoreCaloCut ) {
            // Add fraction m_resNeuhEne of neutron corresponding to EMCal deposit
            tr_ET  +=  m_resNeuhEne * (*p_MC_Iter)->momentum().perp();
            tr_Eta += (*p_MC_Iter)->momentum().pseudoRapidity() * m_resNeuhEne * (*p_MC_Iter)->momentum().perp();
            resphi  = Phi( (*p_MC_Iter)->momentum().phi() - tr_SeedPhi );
            tr_Phi += resphi * m_resNeuhEne * (*p_MC_Iter)->momentum().perp();
          }
        }
      }

      p_log << MSG::DEBUG << "-----> NPi0 = " << tr_NPi0 
            << "\tETPi0 = " <<tr_ETPi0/GeV << endreq; 


      // Check whether candidate has muon or electron tracks...
      // If so, randomly decide whether to flag as muon or 
      // electron veto.
      if( tr_NEleTrack > 0 ){
        float rele_veto = randFlat()->fire();
        if( rele_veto < m_eleVetoEff ) { 
          tr_eleVetoFlag=1;
          p_log << MSG::DEBUG << "Flagged by Electron Veto"<< endreq;
        }
      }
      if( tr_NMuoTrack > 0 ){
        float rmuo_veto = randFlat()->fire();
        if( rmuo_veto < m_muoVetoEff ) { 
          tr_muoVetoFlag=1;
          p_log << MSG::DEBUG << "Flagged by Muon Veto"<< endreq;
        }
      }

      p_log << MSG::DEBUG << tr_NProng << "P Candidate Made."<< endreq;
      p_log << MSG::VERBOSE << "Sum Core Track pT = "<< tr_PTTrack << endreq;
      p_log << MSG::VERBOSE << "Sum Core Neutral pion ET = "<< tr_ETPi0 << endreq;
      p_log << MSG::VERBOSE << "Sum Core ET (including other neutrals) = "<< tr_ET/GeV << endreq;
          
      // Define Bary-Centred eta-phi position...
      tr_Eta = tr_Eta / tr_ET;
      // Get correct -pi<phi<pi value for weighted phi position
      tr_Phi = tr_Phi/tr_ET + tr_SeedPhi;
          
      p_log << MSG::DEBUG << "Candidate eta = "<< tr_Eta <<", phi = "<< tr_Phi << endreq;
          

      //---------------------------------------------------------
      // Check if IsTrueTau
      //---------------------------------------------------------
      std::vector<const HepMC::GenParticle*>::iterator p_partIter;
      for( p_partIter = p_particles.begin(); p_partIter != p_particles.end();
           p_partIter++ )
        {
          //Find Taus in truth record
          if( std::abs( (*p_partIter)->pdg_id() ) == 15 )
            {
            // Find Tau children
              HepMC::GenVertex*  endVertex = (*p_partIter)->end_vertex();
              if(!endVertex) continue;

              HepLorentzVector tauHlv = (*p_partIter)->momentum(); // tau momentum
              HepLorentzVector neuHlv; // neutrino momentum
              HepLorentzVector visHlv; // visible momentum
              
              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 Decaying tau in MC event information
              bool TauDecayVertexInMC = true;         
              for(; thisChild!=endChild; ++thisChild){
                //p_log << MSG::VERBOSE << "     Child pdgid = " << (*thisChild)->pdg_id() <<endreq; 
                if(    abs((*thisChild)->pdg_id()) == 11
                    || abs((*thisChild)->pdg_id()) == 13
                    || abs((*thisChild)->pdg_id()) == 15 )
                  {
                    //If child is a tau then haven't found a hadronic tau decay vertex
                    TauDecayVertexInMC = false;
                    break;
                  }
                // Find Pt of Taus tau-neutrino
                if(  abs((*thisChild)->pdg_id()) == 12 || 
                     abs((*thisChild)->pdg_id()) == 14 || 
                     abs((*thisChild)->pdg_id()) == 16 )
                  {
                    neuHlv += (*thisChild)->momentum();
                  }
              }
              if(TauDecayVertexInMC==false)
                {
                  p_log << MSG::VERBOSE << "Not a Hadronic Tau decay. Ignore Vertex. Move to next Vertex."<<endreq;
                  continue;
                }
              
              visHlv = tauHlv - neuHlv;

              // Check whether matches tau candidate with dR cut
              // dR separation of tau candidate from a true tau
              // should I search for hadronic taus only here??
              Phi    truth_detPhi = tr_Phi - visHlv.phi();
              double truth_detEta = fabs( tr_Eta - visHlv.pseudoRapidity() );
              double truth_detCone = sqrt((truth_detPhi * truth_detPhi) + (truth_detEta * truth_detEta) );
                  
              if ( truth_detCone < m_detRTrueTauCut ){
                p_log << MSG::DEBUG << "Matches tau with visible decay to eta = "
                      << visHlv.pseudoRapidity()
                      <<", phi = "<< visHlv.phi() 
                      << " (dR = "<< truth_detCone << ")" << endreq;
                tr_IsTrueTau=1;
                break;
              }
            }
        }

      //-------------------------------------------------
      // Smear Energy from gaussian distribution
      //-------------------------------------------------         
      float aa, sigma;
      while(1) {
        aa = randGauss()->fire();

        sigma = aa* m_resEflowEne/sqrt( tr_ET * cosh( tr_Eta ) / GeV );
        if( 1.0 +sigma > .0) break;
      }

      p_log << MSG::DEBUG << "After Smearing, ET = "<< tr_ET * (1.0 + sigma) << endreq;

      if ( m_recoETCut > tr_ET * (1.0 + sigma) )
        {
          p_log << MSG::DEBUG << "Tau candidate rejected. (Smeared ET below m_recoETCut = "<<m_recoETCut<<")"<< endreq;
          continue;
        }
      p_log << MSG::DEBUG << "Recording Tau Candidate in Collection"<< endreq;
      //-------------------------------------------------
      // Record new object in collection
      //-------------------------------------------------
      Tau *tr_object = new Tau();
          
      if( tr_PTTrack != 0 )
        {
          tr_EtaTrack /= tr_PTTrack;
          tr_PhiTrack = Phi(tr_PhiTrack/tr_PTTrack + tr_SeedPhi );
        }
          
      if( tr_ETPi0 != 0 )
        {
          tr_EtaPi0 /= tr_ETPi0;
          tr_PhiPi0 = Phi( tr_PhiPi0/tr_ETPi0 + tr_SeedPhi );
        }
          
      // calculate mass of tau candidate
      float pxTrack =  tr_PTTrack * cos( tr_PhiTrack);
      float pxPi0  =  tr_ETPi0  * cos( tr_PhiPi0);
          
      float pyTrack =  tr_PTTrack * sin( tr_PhiTrack);
      float pyPi0  =  tr_ETPi0  * sin( tr_PhiPi0);
          
      float pzTrack =  tr_PTTrack * sinh( tr_EtaTrack);
      float pzPi0  =  tr_ETPi0  * sinh( tr_EtaPi0);
          
      float eeTrack =  tr_PTTrack * cosh( tr_EtaTrack);
      float eePi0  =  tr_ETPi0  * cosh( tr_EtaPi0);
          
      float mass = ( eePi0 + eeTrack )* ( eePi0 + eeTrack )
        - ( pxTrack + pxPi0 )*( pxTrack + pxPi0 )
        - ( pyTrack + pyPi0 )*( pyTrack + pyPi0 )
        - ( pzTrack + pzPi0 )*( pzTrack + pzPi0 );
          
      if( mass > 0 ) { mass = sqrt (mass); }
      if( mass < 0 ) { mass = 0.; }
          
      tr_object->setE( tr_ET * cosh(  tr_Eta )* (1.0 + sigma) );
      tr_object->setM( mass * (1.0 + sigma)  );
      tr_object->setEta( tr_Eta );
      tr_object->setPhi( tr_Phi );
      tr_object->setNTrack( tr_NProng );         
      tr_object->setSumQTrack( tr_SumQ );         
      tr_object->setIsTrueTau( tr_IsTrueTau );

      tr_object->setEleVeto(tr_eleVetoFlag);
      tr_object->setMuoVeto(tr_muoVetoFlag);

      container->push_back( tr_object );
      
    } // end of itr loop 
    
    // Clear Track vectors before next event...
    TrackQualityVec.clear();
    
    // Store the taus (if any) in StoreGate
    stat = m_tesIO->store( container, m_outputLocation );
    std::string mess = stat ? "Store to TES success" : "Store to TES failure";    
    p_log << MSG::DEBUG << mess << endreq;
    
    //---------------------------------------------------------
        
    return StatusCode::SUCCESS;
  }

  StatusCode TauMaker::finalize()
  {

    MsgStream p_log( messageService(), name() ) ;
    p_log << MSG::DEBUG << "TauMaker Finalize "<< endreq;
        
    //   Note from CLHEP/Random/RandFlat instructions:
    //   "If the engine is passed by pointer then the corresponding engine
    //   object will be deleted by the RandFlat destructor.If the engine is
    //   passed by reference then the corresponding engine object will not be
    //   deleted by the RandFlat destructor." 
    //   So don't delete the random engine after m_randflat.

    //still crashing - maybe engine is being deleted when randGauss is deleted,
    // then randFlat is trying to delete the same engine again... can fix by using
    // two different HepRandomEngines. Seems a bit wasteful. 

    //Pass engine to distribution classes by reference rather than pointer. This way the
    //Engine is not deleted when the distribution is deleted

    if (m_randGauss)    delete m_randGauss;
    if (m_randFlat)     delete m_randFlat;
    if (m_randEngine)   delete m_randEngine;


    return StatusCode::SUCCESS;
  }

  
  
  double TauMaker::GetTrackEfficiency( const Track * atlTrack ) {
    // Get efficiency for a given track to pass "good quality"
    // reconstruction cuts. Interpolate as function of pT
    
    double trackEff = 0.0;

    double trackPT = atlTrack->pT();
    double trackID = abs( atlTrack->pdg_id() );

    //Check trackPT is greater than lowest interpolation point
    if( trackPT > m_trackPTPoint.at(0) 
        &&  trackPT < m_trackPTPoint.at( m_trackPTPoint.size()-1 ) ) {
      
      for( unsigned int i=1; i<m_trackPTPoint.size(); i++ ){
        if( m_trackPTPoint.at(i) > trackPT ) {

          double lowPoint;
          double highPoint;

          if(trackID==11){
            lowPoint= m_trackEleEff.at(i-1);
            highPoint=m_trackEleEff.at(i);        
          }else if(trackID==13){
            lowPoint= m_trackMuoEff.at(i-1);
            highPoint=m_trackMuoEff.at(i);        
          }else{
            lowPoint=m_trackHadEff.at(i-1);
            highPoint=m_trackHadEff.at(i);        
          }
          double m = ( lowPoint - highPoint ) / (m_trackPTPoint.at(i-1) - m_trackPTPoint.at(i));
          trackEff = m*trackPT + ( highPoint - m* m_trackPTPoint.at(i) );
          break;
        }
      }
    } else if ( trackPT > m_trackPTPoint.at( m_trackPTPoint.size()-1 ) ) {
      if(trackID==11){
        trackEff = m_trackEleEff.at( m_trackPTPoint.size()-1 );   
      }else if(trackID==13){
        trackEff = m_trackMuoEff.at( m_trackPTPoint.size()-1 );   
      }else{
        trackEff = m_trackHadEff.at( m_trackPTPoint.size()-1 );   
      }
    }
    
    return trackEff;
  }
  
}

---