TauTagger.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------
//
// file:     TauTagger.cxx
//
//-----------------------------------------------------------------------

#include <fstream>
#include <iostream>
#include <vector>
#include <string>
#include <deque>

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

#include "AtlfastEvent/Interpolator.h"

#include "CLHEP/Random/JamesRandom.h"
#include "CLHEP/Random/RandFlat.h"
#include "CLHEP/Random/RandomEngine.h"

#include "CLHEP/Units/SystemOfUnits.h"

#include "PathResolver/PathResolver.h"

namespace Atlfast
{
  
  //------------------------------------------------------------------
  // Constructor
  //------------------------------------------------------------------
  TauTagger::TauTagger( const std::string& name, ISvcLocator* pSvcLocator ) 
    : Algorithm( name, pSvcLocator ),
      m_randEngine(0),
      m_randFlat(0),
      m_tesIO(0),
      m_tau1P_eff_interpolator(0),
      m_tau1P_rej_interpolator(0),
      m_tau2P_eff_interpolator(0),
      m_tau2P_rej_interpolator(0),
      m_tau3P_eff_interpolator(0),
      m_tau3P_rej_interpolator(0)

  {
    //Default private parameters
    
    m_interpolateTauIdent =  true;
    m_tauLocation         =  "/Event/AtlfastTaus";
    m_tau1P_EffFile       =  "atlfastDatafiles/tau1P_EffFile.txt";
    m_tau1P_RejFile       =  "atlfastDatafiles/tau1P_RejFile.txt";
    m_tau2P_EffFile       =  "atlfastDatafiles/tau2P_EffFile.txt";
    m_tau2P_RejFile       =  "atlfastDatafiles/tau2P_RejFile.txt";
    m_tau3P_EffFile       =  "atlfastDatafiles/tau3P_EffFile.txt";
    m_tau3P_RejFile       =  "atlfastDatafiles/tau3P_RejFile.txt";
    
    //Declare joboption properties here...
    declareProperty( "interpolateTauIdent", m_interpolateTauIdent );
    declareProperty( "TauLocation",         m_tauLocation );
    declareProperty( "tau1P_EffFile",       m_tau1P_EffFile );
    declareProperty( "tau1P_RejFile",       m_tau1P_RejFile );
    declareProperty( "tau2P_EffFile",       m_tau2P_EffFile );
    declareProperty( "tau2P_RejFile",       m_tau2P_RejFile );
    declareProperty( "tau3P_EffFile",       m_tau3P_EffFile );
    declareProperty( "tau3P_RejFile",       m_tau3P_RejFile );
  }
  
  TauTagger::~TauTagger() {
    
    if (m_tesIO) {
      delete m_tesIO;
    }
    
  }
  
  //----------------------------------------------------------------
  // Tool initialization
  //----------------------------------------------------------------
  StatusCode TauTagger::initialize()
  {
    MsgStream p_log( msgSvc(), name() );
    
    p_log << MSG::DEBUG << "Initialising TauTagger" << endreq;
    
    GlobalEventData* ged = GlobalEventData::Instance();
    int randSeed = ged->randSeed() ;
    m_tesIO      = new TesIO(ged->mcLocation(), ged->justHardScatter());

    m_randEngine = new HepJamesRandom( randSeed );
    m_randFlat   = new RandFlat( m_randEngine);


    // Create interpolator for tau efficiencies   
    makeInterpolator(m_tau1P_eff_interpolator, m_tau1P_EffFile , true);
    makeInterpolator(m_tau1P_rej_interpolator, m_tau1P_RejFile , true);
    makeInterpolator(m_tau2P_eff_interpolator, m_tau2P_EffFile , true);
    makeInterpolator(m_tau2P_rej_interpolator, m_tau2P_RejFile , true);
    makeInterpolator(m_tau3P_eff_interpolator, m_tau3P_EffFile , true);
    makeInterpolator(m_tau3P_rej_interpolator, m_tau3P_RejFile , true);
    
    return StatusCode::SUCCESS;
  }

  //-----------------------------------------------------------------
  // Tool execution
  //-----------------------------------------------------------------
  StatusCode TauTagger::execute()
  {
    MsgStream p_log( msgSvc(), name() );
        
    p_log << MSG::DEBUG << "Executing TauTagger" << endreq;

    const TauContainer* container(0);
    if(!m_tesIO->getDH(container, m_tauLocation)){
      p_log << MSG::ERROR 
            << "Could not get DataHandle for " << m_tauLocation 
            << endreq ;
      return StatusCode::SUCCESS ;
    }
    
    TauContainer::const_iterator p_iterF = container->begin();
    for( ; p_iterF != container->end(); p_iterF++ )
      {
                
        //-------------------------------------------------
        // Make Identification Decision
        //-------------------------------------------------

        //long tr_NTrack = (*p_iterF)->nTrack();
        //int tr_IsTrueTau = (*p_iterF)->isTrueTau();
        int tr_Discriminant = 0;
        
        double idEff=GetIDEff(p_iterF);
        
        float rid = randFlat()->fire();
        if(rid<idEff) tr_Discriminant = 1;
        p_log << MSG::DEBUG << "Discriminant Tag = "
              << tr_Discriminant <<" ( rand seed = "<< rid 
              << " id Eff. = "<< idEff << " )" << endreq;       
        (*p_iterF)->setDiscriminant(tr_Discriminant);
      }
    
    return StatusCode::SUCCESS;
  }


  StatusCode TauTagger::finalize()
  {
    
    MsgStream p_log( msgSvc(), name() );        
    p_log << MSG::DEBUG << "TauTagger 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.

    if (m_randFlat)      delete m_randFlat;

    if (m_tau1P_eff_interpolator)      delete m_tau1P_eff_interpolator;
    if (m_tau1P_rej_interpolator)      delete m_tau1P_rej_interpolator;
    if (m_tau2P_eff_interpolator)      delete m_tau2P_eff_interpolator;
    if (m_tau2P_rej_interpolator)      delete m_tau2P_rej_interpolator;
    if (m_tau3P_eff_interpolator)      delete m_tau3P_eff_interpolator;
    if (m_tau3P_rej_interpolator)      delete m_tau3P_rej_interpolator;

    return StatusCode::SUCCESS;
  }



  double TauTagger::GetIDEff(TauContainer::const_iterator p_iterF){

    // Find Identification efficiency for tau candidate to pass
    // discriminant.

    MsgStream p_log( msgSvc(), name() );

    double IDEff = 0.0;

    // Input to interpolators
    double pt_in_GeV = (*p_iterF)->pt()/GeV;
    double etatau = (*p_iterF)->eta();
    int nTrack = (int) (*p_iterF)->nTrack();
    deque<double> input_values;
    input_values.push_back( pt_in_GeV );
    input_values.push_back( fabs(etatau) );

    if (nTrack==1){
      if( (*p_iterF)->isTrueTau()==1 ){
        IDEff = m_interpolateTauIdent ?
          m_tau1P_eff_interpolator->interpolate(input_values) :
          m_tau1P_eff_interpolator->getNearestValue(input_values);
      } else if ( (*p_iterF)->isTrueTau()==0 ){
        IDEff = m_interpolateTauIdent ?
          m_tau1P_rej_interpolator->interpolate(input_values) :
          m_tau1P_rej_interpolator->getNearestValue(input_values);
      }
    }
    if (nTrack==2){
      if( (*p_iterF)->isTrueTau()==1 ){
        IDEff = m_interpolateTauIdent ?
          m_tau2P_eff_interpolator->interpolate(input_values) :
          m_tau2P_eff_interpolator->getNearestValue(input_values);
      } else if ( (*p_iterF)->isTrueTau()==0 ){
        IDEff = m_interpolateTauIdent ?
          m_tau2P_rej_interpolator->interpolate(input_values) :
          m_tau2P_rej_interpolator->getNearestValue(input_values);
      }
    }
    if (nTrack>=3){
      if( (*p_iterF)->isTrueTau()==1 ){
        IDEff = m_interpolateTauIdent ?
          m_tau3P_eff_interpolator->interpolate(input_values) :
          m_tau3P_eff_interpolator->getNearestValue(input_values);
      } else if ( (*p_iterF)->isTrueTau()==0 ){
        IDEff = m_interpolateTauIdent ?
          m_tau3P_rej_interpolator->interpolate(input_values) :
          m_tau3P_rej_interpolator->getNearestValue(input_values);
      }
    }

    p_log << MSG::DEBUG << "Tau"<<nTrack<<"P true="<<(*p_iterF)->isTrueTau()
          <<";  pT= "<<pt_in_GeV << " GeV; eta= "<<etatau 
          <<";\tideff= "<<IDEff <<endreq;

    return IDEff;
  }


  void TauTagger::makeInterpolator(Interpolator* &intptr, std::string intname, bool contbounds){
    MsgStream p_log( msgSvc(), name() );
    std::string fn = PathResolver::find_file(intname, "DATAPATH");
    intptr = new Interpolator(fn);
    p_log << MSG::DEBUG << "Made Interpolator from input file " << fn << endreq;
    if (contbounds) intptr->setContinuousBoundaries();
    return;
  }

}

---