Atlfast::TauMaker Class Reference

Tool to create reference TauFastObject (for physics sample). More...

#include <TauMaker.h>

Collaboration diagram for Atlfast::TauMaker:

[legend]List of all members.

Public Member Functions
	TauMaker (const std::string &name, ISvcLocator *pSvcLocator)
	~TauMaker ()
StatusCode	initialize ()
StatusCode	execute ()
StatusCode	finalize ()
RandGauss *	randGauss ()
RandFlat *	randFlat ()
double	GetTrackEfficiency (const Track *)
Private Attributes
TesIO *	m_tesIO
HepRandomEngine *	m_randEngine
RandGauss *	m_randGauss
RandFlat *	m_randFlat
std::vector< double >	m_trackPTPoint
std::vector< double >	m_trackHadEff
std::vector< double >	m_trackMuoEff
std::vector< double >	m_trackEleEff
std::string	m_TrackCollection
	name of the atlfast Track container
std::string	m_outputLocation
	name of the output tau collection
double	m_recoEtaCut
	rapidity coverage for reconstruction
double	m_pTLeadTrackCut
	transverse momenta threshold for leading track
double	m_pTOtherTrackCut
	transverse momenta threshold for other tracks
double	m_maxNumTracks
	multiplicity of nearby tracks (qualified)
double	m_detRCoreTrackCut
	half-angle opening for nearby tracks
double	m_detRCoreCaloCut
	half-angle opening for calorimetric signal region
double	m_detRRecoSeparation
	half-angle opening for seed track from other reconstructed candidate
double	m_detRTrueTauCut
	half-angle opening truth matching
double	m_resEflowEne
	assumed resolution for energy smearing
double	m_recoETCut
	Cut on reconstructed ET.
double	m_resNeuhEne
	assumed EM response for neutrons
double	m_eleVetoEff
	assumed efficiency for vetoing candidates from electrons
double	m_muoVetoEff
	assumed efficiency for vetoing candidates from muons
double	m_resTrack
	assumed EM response for charged pion tracks
std::string	m_TrackEffFile
	Track Efficiency Parametrisation - Input File.

---

Detailed Description

Tool to create reference TauFastObject (for physics sample).

Loops over all particles in collection to find true tau, selects hadronic decays calculates ET, (eta, phi) of the visible part, multiplicity of charged (pi+-) and neutral (pi0) components, bary-center of charged and neutral system

Definition at line 29 of file TauMaker.h.

---

Constructor & Destructor Documentation

Atlfast::TauMaker::TauMaker	(	const std::string &	name,
		ISvcLocator *	pSvcLocator
	)

Constructor

Parameters:

	type	type of tool
	name	name of tool
	parent	pointer to the IInterface of the parent component

Definition at line 38 of file TauMaker.cxx.

    : 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       );
  }

Atlfast::TauMaker::~TauMaker

(

)

Definition at line 84 of file TauMaker.cxx.

                      {
    
    if (m_tesIO) {
      delete m_tesIO;
    }
  }

---

Member Function Documentation

StatusCode Atlfast::TauMaker::initialize

(

)

Algorithm initialization

Returns:

Status code

Definition at line 94 of file TauMaker.cxx.

  {
    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;
  }

StatusCode Atlfast::TauMaker::execute

(

)

Algorithm execution

Returns:

Status code

Definition at line 170 of file TauMaker.cxx.

  {

    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 Atlfast::TauMaker::finalize

(

)

Algorithm finalization

Returns:

Status code

Definition at line 775 of file TauMaker.cxx.

  {

    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;
  }

RandGauss* Atlfast::TauMaker::randGauss

(

)

[inline]

Definition at line 65 of file TauMaker.h.

{ return m_randGauss;}

RandFlat* Atlfast::TauMaker::randFlat

(

)

[inline]

Definition at line 66 of file TauMaker.h.

{ return m_randFlat;}

double Atlfast::TauMaker::GetTrackEfficiency

(

const Track *

)

Definition at line 805 of file TauMaker.cxx.

                                                              {
    // 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;
  }

---

Member Data Documentation

TesIO* Atlfast::TauMaker::m_tesIO [private]

Definition at line 72 of file TauMaker.h.

HepRandomEngine* Atlfast::TauMaker::m_randEngine [private]

Definition at line 74 of file TauMaker.h.

RandGauss* Atlfast::TauMaker::m_randGauss [private]

Definition at line 75 of file TauMaker.h.

RandFlat* Atlfast::TauMaker::m_randFlat [private]

Definition at line 76 of file TauMaker.h.

std::vector<double> Atlfast::TauMaker::m_trackPTPoint [private]

Definition at line 80 of file TauMaker.h.

std::vector<double> Atlfast::TauMaker::m_trackHadEff [private]

Definition at line 81 of file TauMaker.h.

std::vector<double> Atlfast::TauMaker::m_trackMuoEff [private]

Definition at line 82 of file TauMaker.h.

std::vector<double> Atlfast::TauMaker::m_trackEleEff [private]

Definition at line 83 of file TauMaker.h.

std::string Atlfast::TauMaker::m_TrackCollection [private]

name of the atlfast Track container

Definition at line 88 of file TauMaker.h.

std::string Atlfast::TauMaker::m_outputLocation [private]

name of the output tau collection

Definition at line 90 of file TauMaker.h.

double Atlfast::TauMaker::m_recoEtaCut [private]

rapidity coverage for reconstruction

Definition at line 92 of file TauMaker.h.

double Atlfast::TauMaker::m_pTLeadTrackCut [private]

transverse momenta threshold for leading track

Definition at line 94 of file TauMaker.h.

double Atlfast::TauMaker::m_pTOtherTrackCut [private]

transverse momenta threshold for other tracks

Definition at line 96 of file TauMaker.h.

double Atlfast::TauMaker::m_maxNumTracks [private]

multiplicity of nearby tracks (qualified)

Definition at line 98 of file TauMaker.h.

double Atlfast::TauMaker::m_detRCoreTrackCut [private]

half-angle opening for nearby tracks

Definition at line 100 of file TauMaker.h.

double Atlfast::TauMaker::m_detRCoreCaloCut [private]

half-angle opening for calorimetric signal region

Definition at line 102 of file TauMaker.h.

double Atlfast::TauMaker::m_detRRecoSeparation [private]

half-angle opening for seed track from other reconstructed candidate

Definition at line 104 of file TauMaker.h.

double Atlfast::TauMaker::m_detRTrueTauCut [private]

half-angle opening truth matching

Definition at line 106 of file TauMaker.h.

double Atlfast::TauMaker::m_resEflowEne [private]

assumed resolution for energy smearing

Definition at line 108 of file TauMaker.h.

double Atlfast::TauMaker::m_recoETCut [private]

Cut on reconstructed ET.

Definition at line 110 of file TauMaker.h.

double Atlfast::TauMaker::m_resNeuhEne [private]

assumed EM response for neutrons

Definition at line 112 of file TauMaker.h.

double Atlfast::TauMaker::m_eleVetoEff [private]

assumed efficiency for vetoing candidates from electrons

Definition at line 115 of file TauMaker.h.

double Atlfast::TauMaker::m_muoVetoEff [private]

assumed efficiency for vetoing candidates from muons

Definition at line 117 of file TauMaker.h.

double Atlfast::TauMaker::m_resTrack [private]

assumed EM response for charged pion tracks

Definition at line 120 of file TauMaker.h.

std::string Atlfast::TauMaker::m_TrackEffFile [private]

Track Efficiency Parametrisation - Input File.

Definition at line 123 of file TauMaker.h.

---

The documentation for this class was generated from the following files:

• /afs/cern.ch/user/s/sdean/project_builds/AtlasOffline-rel_1/Simulation/Atlfast/AtlfastAlgs/AtlfastAlgs/TauMaker.h
• /afs/cern.ch/user/s/sdean/project_builds/AtlasOffline-rel_1/Simulation/Atlfast/AtlfastAlgs/src/TauMaker.cxx

---