AtlfastB.cxx

Go to the documentation of this file.

// ================================================
// AtlfastB class Implementation
// ================================================
//
// THIS TEXT TO BE REPLACED BY ATLAS STANDARD FORMAT
//
// Namespace ATLF::
//

#include <fstream>
#include <cmath>
#include <deque>
#include <utility>

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

#include "AtlfastEvent/Jet.h"
#include "AtlfastEvent/IKinematic.h"
#include "AtlfastEvent/Interpolator.h"

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

// Gaudi includes
#include "GaudiKernel/DataSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"

//Other Packages used by this class
// Other Packages used by this class:-

#include "CLHEP/Random/Ranlux64Engine.h"
#include "CLHEP/Random/RandFlat.h"
#include "CLHEP/Units/SystemOfUnits.h"

#include "PathResolver/PathResolver.h"

//--------------------------------
// Constructors and destructors
//--------------------------------

namespace Atlfast {
  
  using std::abs;

  AtlfastB::AtlfastB( const std::string& name, ISvcLocator* pSvcLocator ) 
    : Algorithm( name, pSvcLocator ),
      m_epsilonBjet(0),
      m_beff_interpolator(0),
      m_brejpu_interpolator(0),
      m_brejnpu_interpolator(0),
      m_brejtau_interpolator(0),
      m_brejc_interpolator(0),
      m_tesIO(0),
      m_pRandomEngine(0),
      m_pRandFlatGenerator(0),
      m_Jets(0),
      m_taueff_interpolator(0),
      m_taurej_interpolator(0)
  {

    MsgStream log( messageService(), name ) ;
    log << MSG::DEBUG << "Constructor" << endreq;
    
    //Default private parameters
    
    m_AtlfBJetSwitch      =  true;
    m_AtlfCalSwitch       =  true;
    m_AtlfTauSwitch       =  true;
    m_AtlfTauVetoSwitch   =  false;
    m_AtlfTrigMuoSwitch   =  false;
    m_atlfBNSet           =  10;
    m_epsitau             =  0.5;
    m_indtauveto          =  1;
    m_corrjfile           =  "atlfastDatafiles/AtlfastBjet.dat";
    m_corrcfile           =  "atlfastDatafiles/AtlfastBcjet.dat";
    m_corrtaumom          =  false;
    m_useTDRBParam        =  false;
    m_interpolateBTagging =  false;
    m_correctionFactor    =  1.;

    m_AtlfTau1P3PSwitch   =  false;
    m_epsitau1P           =  0.35;
    m_epsitau3P           =  0.08;
    m_Tau1P3Pcorrfile    =  "atlfastDatafiles/Tau1P3Pcorrfile.txt";

    // Default paths for entities in the TES
    m_inputLocation     = "/Event/AtlfastJets";

    // This is how you declare the paramemters to Gaudi so that
    // they can be over-written via the job options file
    
    
    
    
    declareProperty( "InputLocation",        m_inputLocation ) ;
    
    declareProperty( "AtlfBjetSwitch",       m_AtlfBJetSwitch ) ;
    declareProperty( "AtlfCalSwitch",        m_AtlfCalSwitch ) ;
    declareProperty( "AtlfTauSwitch",        m_AtlfTauSwitch );
    declareProperty( "AtlfTauVetoSwitch",    m_AtlfTauVetoSwitch );
    declareProperty( "AtlfTrigMuoSwitch",    m_AtlfTrigMuoSwitch );
    declareProperty( "AtlfTau1P3PSwitch",    m_AtlfTau1P3PSwitch );

    
    declareProperty( "AtlfBNSet",            m_atlfBNSet ) ;
    declareProperty( "TauEff",               m_epsitau ) ;
    declareProperty( "TauVetoOption",        m_indtauveto ) ;
    declareProperty( "Tau1PEff",             m_epsitau1P ) ;
    declareProperty( "Tau3PEff",             m_epsitau3P ) ;
    
    declareProperty( "JetCorrFile",          m_corrjfile ) ;
    declareProperty( "CJetCorrFile",         m_corrcfile ) ;
    declareProperty( "Tau1P3PCorrFile",      m_Tau1P3Pcorrfile ) ;
    
    declareProperty( "UseTDRBParam",         m_useTDRBParam );
    declareProperty( "InterpolateBTagging",  m_interpolateBTagging );
    declareProperty( "BCorrectionFactor",    m_correctionFactor); 
    // m_correctionFactor less than 1 : degradation of peformances

    declareProperty( "CalibrateTauMomentum", m_corrtaumom );

    log << MSG::DEBUG << "End of constructor" << endreq;

  }
  
  AtlfastB::~AtlfastB() {
    
    MsgStream log( messageService(), name() ) ;
    log << MSG::DEBUG << "Destructor" << endreq;
    if ( m_tesIO )                delete m_tesIO;
    if ( m_pRandomEngine )        delete m_pRandomEngine;
    if ( m_pRandFlatGenerator )   delete m_pRandFlatGenerator;
    if ( m_taueff_interpolator )  delete m_taueff_interpolator;
    if ( m_taurej_interpolator )  delete m_taurej_interpolator;
    if ( m_beff_interpolator )    delete m_beff_interpolator;
    if ( m_brejpu_interpolator )  delete m_brejpu_interpolator;
    if ( m_brejnpu_interpolator ) delete m_brejnpu_interpolator;
    if ( m_brejtau_interpolator ) delete m_brejtau_interpolator;
    if ( m_brejc_interpolator )   delete m_brejc_interpolator;
  }
  
  
  
  //---------------------------------
  // initialise() 
  //---------------------------------
  
  StatusCode AtlfastB::initialize(){
    MsgStream log( messageService(), name() ) ;
    log << MSG::DEBUG << "instantiating an AtlfastB" << endreq;
    

    m_corrjfile         =  PathResolver::find_file (m_corrjfile, "DATAPATH");
    m_corrcfile         =  PathResolver::find_file (m_corrcfile, "DATAPATH");
    m_Tau1P3Pcorrfile   =  PathResolver::find_file (m_Tau1P3Pcorrfile, "DATAPATH");

    log << MSG::DEBUG << " m_corrjfile after PathResolver: " << m_corrjfile <<endreq;
    log << MSG::DEBUG << " m_corrcfile after PathResolver: " << m_corrcfile <<endreq;
    log << MSG::DEBUG << " m_Tau1P3Pcorrfile after PathResolver: " << m_Tau1P3Pcorrfile <<endreq;

    //Had to move this section above the TesIO call
    //Flat random number generator
    //get the Global Event Data using singleton pattern

    GlobalEventData* ged = GlobalEventData::Instance();
    int randSeed = ged->randSeed() ;
    // load the location of the MC in StoreGate
    m_mcLocation       = ged -> mcLocation();

    m_tesIO= new TesIO(m_mcLocation, ged->justHardScatter());
    
    std::string fn; // For filenames

    if (m_useTDRBParam) {
      //open files for jet correction factors
      
      std::ifstream inputjetfile;
      inputjetfile.open(m_corrjfile.c_str());
      
      if(inputjetfile){
        
        log << MSG::INFO 
            << "Pt jet correction factors file " 
            << m_corrjfile
            <<" open."
            <<endreq;
        int nrow;
        int ncolumn;
        inputjetfile>>nrow;
        inputjetfile>>ncolumn;
        if(nrow != 5||ncolumn != 8) {
          log << MSG::ERROR
              <<"no. of input rows,columns: "
              <<nrow<<" "
              <<ncolumn<<" "
              <<"expected 5, 8 "
              <<endreq;  
          return StatusCode::FAILURE ;
        }
        
        for(int i=0; i<nrow;i++){
          for(int j=0; j<ncolumn; j++){
            
            int nset = -1;
            switch (j){
              
            case 0:
              nset=1;
              break;
            case 1:
              nset=2;
              break;
            case 2:
              nset=3;
              break;
            case 3:
              nset=5;
              break;
            case 4:
              nset=11;
              break;
            case 5:
              nset=12;
              break;
            case 6:
              nset=13;
              break;
            case 7:
              nset=14;
              break;
            default:
              log << MSG::ERROR<<"j="<<j<<" Nset= "<<nset
                  <<"  value not correct. Nset must b one of those: 1,2,3,5,11,12,13,14 "<<endreq;
              return StatusCode::FAILURE;
              break;
              
              
            }
            inputjetfile>>m_corrj[i][nset];
            
          }
        }
        
      }else{
        log << MSG::ERROR << "Pt jet correction factors file not found" <<endreq;
        return StatusCode::FAILURE ;
      }
      inputjetfile.close();
      
      std::ifstream inputcjetfile;
      inputcjetfile.open(m_corrcfile.c_str());
      
      if(inputcjetfile){
        
        log << MSG::INFO 
            << "Pt c jet correction factors file " 
            << m_corrcfile
            <<" open."
            <<endreq;
        int nrow;
        int ncolumn;
        inputcjetfile>>nrow;
        inputcjetfile>>ncolumn;
        if(nrow != 5||ncolumn != 8) {
          log << MSG::ERROR
              <<"no. of input rows,columns: "
              <<nrow<<" "
              <<ncolumn<<" "
              <<"expected 5, 8 "
              <<endreq;  
          return StatusCode::FAILURE ;
        }
        for(int i=0; i<nrow;i++){
          for(int j=0; j<ncolumn; j++){
            
            int nset;
            switch (j){
              
            case 0:
              nset=1;
              break;
            case 1:
              nset=2;
              break;
            case 2:
              nset=3;
              break;
            case 3:
              nset=5;
              break;
            case 4:
              nset=11;
              break;
            case 5:
              nset=12;
              break;
            case 6:
              nset=13;
              break;
            case 7:
              nset=14;
              break;
            default:
              log << MSG::ERROR<<"Nset value not correct. Nset must be one of these: 1,2,3,5,11,12,13,14 "<<endreq;
              return StatusCode::FAILURE;
              break;
            }
            inputcjetfile>>m_corrc[i][nset];
          }
        }
        
      }else{
        log << MSG::ERROR 
            << "Pt c jet correction factors file not found" 
            <<endreq;
        return StatusCode::FAILURE ;
      }
      inputcjetfile.close();
    } else {
      int ialgo = m_atlfBNSet; 
      if (ialgo < 1 || (ialgo > 14 && ialgo != 100) ) {
        log << MSG::FATAL << "NSET not correct. Should be between 1 and 14, or 100 for canonical b-tagging" << endreq;
        return StatusCode::FAILURE;
      }
      
      // Make new Interpolators for the b-tagging efficiencies and rejections
      std::string path = "atlfastDatafiles/RejectionFactor";
      std::ostringstream o;
      o << m_atlfBNSet;
      
      makeInterpolator(m_beff_interpolator,"atlfastDatafiles/BEfficiency"+o.str()+".txt",true);
      makeInterpolator(m_brejpu_interpolator,path+"PU"+o.str()+".txt",true);
      makeInterpolator(m_brejnpu_interpolator,path+"NPU"+o.str()+".txt",true);
      makeInterpolator(m_brejtau_interpolator,path+"Tau"+o.str()+".txt",true);
      makeInterpolator(m_brejc_interpolator,path+"C"+o.str()+".txt",true);


      // The efficiency is defined by m_atlfBNSet (for m_atlfBNSet <= 14)
      std::string algName = "IP2D"; 
      if (ialgo == 100) {
        algName = "Canonical";
        m_epsilonBjet = 0.6;
      } else {
        if (ialgo > 7){
          algName = "SV1+IP3D";
          ialgo -= 7;
        }
        m_epsilonBjet = 0.5 + (ialgo - 1)*0.05;
      }
      
      log << MSG::INFO << "Running with " << algName << " at eps_b = " << m_epsilonBjet << endreq;
      
    }
    
    // Read Tau1P3P energy rescaling information from file
    std::ifstream inputfileTau1P3P;
    inputfileTau1P3P.open(m_Tau1P3Pcorrfile.c_str());
    
    if(inputfileTau1P3P){

      int nrow;
      int ncolumn;
      inputfileTau1P3P>>nrow;
      inputfileTau1P3P>>ncolumn;

      // Read corr factor bins from file
      for(int i=0; i<ncolumn; i++){inputfileTau1P3P>>m_corr1P3P[i];}

      // Read corr factor probabilities for Tau1P in 6 PT bins
      for(int j=0; j<nrow; j++){
        for(int k=0; k<ncolumn; k++){
          inputfileTau1P3P>>m_corr_prob1P[j][k];
        }
      }    
      // Read corr factor probabilities for Tau3P in 6 PT bins
      for(int l=0; l<nrow; l++){
        for(int m=0; m<ncolumn; m++){
          inputfileTau1P3P>>m_corr_prob3P[l][m];
        }
      }
    }else{
      log << MSG::ERROR 
          << "Tau1P3P corr factor file not found" 
          <<endreq;
      return StatusCode::FAILURE ;
    }

    // Check Tau1P3P efficiencies are within range
    if( m_epsitau1P < 0.20){m_epsitau1P = 0.00;
    log << MSG::INFO<< "Tau1PEff below allowed range (<0.10). Set Tau1PEff = 0.0 "<<endreq;}
    if( m_epsitau1P > 0.40){m_epsitau1P = 0.40;
    log << MSG::INFO<< "Tau1PEff above allowed range (>0.40). Set Tau1PEff = 0.40"<<endreq;}
    if( m_epsitau3P < 0.05){m_epsitau3P = 0.00;
    log << MSG::INFO<< "Tau3PEff below allowed range (<0.05). Set Tau3PEff = 0.0 "<<endreq;}
    if( m_epsitau3P > 0.10){m_epsitau3P = 0.10;
    log << MSG::INFO<< "Tau3PEff above allowed range (>0.10). Set Tau3PEff = 0.10"<<endreq;}

    
    m_pRandomEngine = new Ranlux64Engine(randSeed);
    m_pRandFlatGenerator=new RandFlat(*m_pRandomEngine);
    
    // Create interpolator for tau jet efficiencies
    
    makeInterpolator(m_taueff_interpolator,"atlfastDatafiles/TauEfficiencies.txt",false);
    makeInterpolator(m_taurej_interpolator,"atlfastDatafiles/TauRejections.txt",false);

    HeaderPrinter hp("AtlfastB:", log);
    
    hp.add("TES Locations:              ");
    hp.add(" Jets from                     ", m_inputLocation); 
    hp.add("AtlfBje switch:                ", m_AtlfBJetSwitch);
    hp.add("AtlfBje NSet:                  ", m_atlfBNSet);
    hp.add("AtlfCal switch:                ", m_AtlfCalSwitch);
    hp.add("AtlfTau switch:                ", m_AtlfTauSwitch);    
    hp.add("AtlfTauVeto switch:            ", m_AtlfTauVetoSwitch);
    hp.add("AtlfTrigMuo switch:            ", m_AtlfTrigMuoSwitch); 
    hp.add("AtlfTau1P3P switch:            ", m_AtlfTau1P3PSwitch);
    hp.add("correction file 1:             ", m_corrjfile); 
    hp.add("correction file 2:             ", m_corrcfile); 
    hp.print();
    
    return StatusCode::SUCCESS ;
  }
  
  //---------------------------------
  // finalise() 
  //---------------------------------
  
  StatusCode AtlfastB::finalize(){
    
    MsgStream log( messageService(), name() ) ;
    log << MSG::INFO << "finalizing" << endreq;    
    return StatusCode::SUCCESS ;
  }
  
  
  //----------------------------------------------
  // execute() method called once per event
  //----------------------------------------------
  
  StatusCode AtlfastB::execute(){
    
    StatusCode sc;
    MsgStream log( messageService(), name() ) ;  
    
    log << MSG::DEBUG<<"In execute"<<endreq;

    if(!m_tesIO->getDH(m_Jets, m_inputLocation)){
      log << MSG::ERROR
          << "Couldn't retrieve JetCollection" << m_inputLocation 
          << endreq ;
      return StatusCode::FAILURE;
    }
    
    if(m_AtlfBJetSwitch){
      //Apply efficiency for b-jet identification
      

      
      sc = m_useTDRBParam ? atlfBje_TDR() : atlfBje();
      
      if(sc.isFailure()){
        log << MSG::ERROR<<"Error in atlfBje"<<endreq;
        return StatusCode::FAILURE;
      }
 
    }

    if(m_AtlfTauSwitch){
      
      //Apply efficiency for tau-jet identification
      
      sc=atlfTau(m_epsitau);
      
      if(sc.isFailure()){
        log << MSG::ERROR<<"Error in atlfTau"<<endreq;
        return StatusCode::FAILURE;
      }

    }
    
    if(m_AtlfTau1P3PSwitch){

      //Apply Tau1P3P efficiency for tau-jet identification

      sc=atlfTau1P3P();

      if(sc.isFailure()){
        log << MSG::ERROR<<"Error in atlfTau1P3P"<<endreq;
        return StatusCode::FAILURE;
      }

    }

    if(m_AtlfTauVetoSwitch){
      
      //Apply veto for tau-jet
      
      sc=atlfTauVeto(m_indtauveto);
      
      if(sc.isFailure()){
        log << MSG::ERROR<<"Error in atlfTauVeto"<<endreq;
        return StatusCode::FAILURE;
      }

    }
    
 
    if(m_AtlfTrigMuoSwitch){
      
      //Apply efficiency for muon trigger. Not implemented yet!
      
      sc=atlfTrigMuo();
      
      if(sc.isFailure()){
        log << MSG::ERROR<<"Error in atlfTrigMuo"<<endreq;
        return StatusCode::FAILURE;
      }
      
    }   
    
    
    if(m_AtlfCalSwitch){
      
      //Recalibrate jet energies
      
      sc=atlfCal();
      
      if(sc.isFailure()){
        log << MSG::ERROR<<"Error in atlfCal"<<endreq;
        return StatusCode::FAILURE;
      }
      
    }
    
    
    //......................................
    // Fill jets that have been successfully tagged and calibrated in JetCollection 
    //and register them in the transient event store. 
    
    
    /*
    JetCollection* myCalJets=new JetCollection;
    
    for(m_it=m_Jets.begin();m_it<m_Jets.end();++m_it){
      myCalJets->push_back(*m_it);
    }

    */
    
    //JetCollection::const_iterator jetcolit = m_Jets->begin();
    m_it = m_Jets->begin();
    

    for( ; m_it != m_Jets->end(); ++m_it ){

      log << MSG::DEBUG
          <<"original jets: PDGID= "<<(*m_it)->pdg_id()
          <<" B tag= "   <<(*m_it)->isBTagged()
          <<" B tag corr factor = " << (*m_it)->bTagCorrFactor()
          <<" Tau tag= " <<(*m_it)->isTauTagged()
          <<" Tau tag corr factor = " << (*m_it)->tauTagCorrFactor()
          <<" pT = "     << (*m_it)->pT() << endreq;
      /*
      log << MSG::DEBUG
          <<"Jets After AtlfastB: PDGID= "<<(*jetcolit)->pdg_id()
          <<" B tag= "    <<(*jetcolit)->isBTagged()
          <<" B tag corr factor = " << (*m_it)->bTagCorrFactor()
          <<" Tau tag= "  <<(*jetcolit)->isTauTagged()
          <<" Tau tag corr factor = " << (*m_it)->tauTagCorrFactor()
          <<" pT = " << (*jetcolit)->pT() << endreq;
      */
      //++m_it;
    }
    

    
    TesIoStat stat = m_tesIO->lock( m_Jets ) ;
    if(!stat){
      log<<MSG::ERROR<<"Could not lock jet collection"<<endreq;
      return stat;
    }
    
    
    return StatusCode::SUCCESS; 
  }
  
  //-------------------------------
  // helper methods implementation
  //-------------------------------
  
  StatusCode AtlfastB::atlfBje_TDR(){
    
    MsgStream log( messageService(), name() ) ; 
    

    switch(m_atlfBNSet){
      
    case 1:
      m_epsib=0.5;
      m_epsic=1./10.9;
      m_epsij=1./231.;
      break;
      
    case 2:
      m_epsib=0.60;
      m_epsic=1./6.7;
      m_epsij=1./93.;
      break;
      
    case 3:
      m_epsib=0.70;
      m_epsic=1./4.3;
      m_epsij=1./34.1;
      break;
      
    case 5:
      m_epsib=0.60;
      m_epsic=1./10.;
      m_epsij=1./100.;
      break;
      
    case 11:
      m_epsib=0.33;
      m_epsic=1./22.9;
      m_epsij=1./1381.;
      break;
      
    case 12:
      m_epsib=0.43;
      m_epsic=1./10.8;
      m_epsij=1./219.;
      break;
      
    case 13:
      m_epsib=0.53;
      m_epsic=1./6.7;
      m_epsij=1./91.;
      break;
      
    case 14:
      m_epsib=0.624;
      m_epsic=1./6.7;
      m_epsij=1./91.;
      break;
      
    default:
      log << MSG::ERROR<<"No b jet tagging efficiencies applied"<<endreq;
      return StatusCode::FAILURE;
      break;
    }
    
    
    //Apply randomized B tagging
    
    log<<MSG::DEBUG<<"Applying randomized B tagging...."<<endreq;
    
    for(m_it=m_Jets->begin();m_it<m_Jets->end();++m_it){
      
      double randnum;
      double corpt;    
      double corr;
      
      double ptjet=(*m_it)->pT();
      double etajet=(*m_it)->eta();
      int pdgid=(*m_it)->pdg_id();
      
      //apply b tagging as for recalibrated jets
      
      
      corr=fitcoreb(ptjet);

      ptjet=corr*ptjet;

      randnum=m_pRandFlatGenerator->shoot(m_pRandomEngine);
      
      corpt=0.;
      
      if(abs(pdgid)==5){

        if((abs(etajet)<=2.5)&&(randnum<m_epsib)){
          (*m_it)->setBTagged("TDR");
        }
        
      }else if(abs(pdgid)==4){
        
        if(ptjet<=30.*GeV) corpt=m_corrc[0][m_atlfBNSet];
        if((ptjet>30.*GeV)&&(ptjet<= 45.*GeV)) corpt=m_corrc[1][m_atlfBNSet];
        if((ptjet>45.*GeV)&&(ptjet<= 60.*GeV)) corpt=m_corrc[2][m_atlfBNSet];   
        if((ptjet>60.*GeV)&&(ptjet<=100.*GeV)) corpt=m_corrc[3][m_atlfBNSet];
        if(ptjet>100.*GeV) corpt=m_corrc[4][m_atlfBNSet];

        if((abs(etajet)<=2.5)&&(randnum<(m_epsic/corpt))) {
          (*m_it)->setBTagged("TDR");
        }
        
      }else if(((abs(pdgid)!=4)&&(abs(pdgid)!=5))||abs(pdgid)==15){
        
        if(ptjet<=30.*GeV) corpt=m_corrj[0][m_atlfBNSet];
        if((ptjet>30.*GeV)&&(ptjet<=45.*GeV)) corpt=m_corrj[1][m_atlfBNSet];
        if((ptjet>45.*GeV)&&(ptjet<=60.*GeV)) corpt=m_corrj[2][m_atlfBNSet];    
        if((ptjet>60.*GeV)&&(ptjet<=100.*GeV)) corpt=m_corrj[3][m_atlfBNSet];
        if(ptjet>100.*GeV) corpt=m_corrj[4][m_atlfBNSet];

        if((abs(etajet)<=2.5)&&(randnum<(m_epsij/corpt))) {
          (*m_it)->setBTagged("TDR");
        }

      }
    } 

    return StatusCode::SUCCESS;

  }
  
  StatusCode AtlfastB::atlfBje() {
    
    MsgStream mlog(messageService(), name());
    //
    mlog << MSG::DEBUG << "Applying randomized B tagging, new parametrization" <<endreq;

    int debug = 0;
    if (mlog.level() < MSG::DEBUG) debug = 1;
    
    for(m_it=m_Jets->begin(); m_it<m_Jets->end(); ++m_it){
      
      double randnum = m_pRandFlatGenerator->shoot(m_pRandomEngine);   
      double ptjet   = (*m_it)->pT();
      double etajet  = (*m_it)->eta();
      int pdgid      = (*m_it)->pdg_id();

      // Input to interpolators
      deque<double> input_values;
      input_values.push_back(ptjet/1.e3);
      input_values.push_back(fabs(etajet));

      if ( fabs(etajet) <= 2.5 ) {
        if (abs(pdgid) == 5) {
          // useless for NSET <= 14 (flat fixed efficiency) and 100 (canonical)
          // useful for NSET > 14 (!=100) (fixed cut)
          double epsilonBjet = m_interpolateBTagging ?
            m_beff_interpolator->interpolate(input_values) :
            m_beff_interpolator->getNearestValue(input_values);
          mlog << MSG::DEBUG << " epsb( "<<ptjet/1.e3<<" , "<<etajet<<" ) = "<< epsilonBjet << endreq;
          if ( randnum < epsilonBjet )
            (*m_it)->setBTagged("SV1+IP3D");
        } else if (abs(pdgid) == 4) {
          double RejC = m_interpolateBTagging ?
            m_brejc_interpolator->interpolate(input_values) :
            m_brejc_interpolator->getNearestValue(input_values);
          mlog << MSG::DEBUG << " Rc( "<<ptjet/1.e3<<" , "<<etajet<<" ) = "<< RejC << endreq;
          if ( randnum < 1./RejC )
            (*m_it)->setBTagged("SV1+IP3D");
        } else if (abs(pdgid) == 15) {
          double RejTau = m_interpolateBTagging ?
            m_brejtau_interpolator->interpolate(input_values) :
            m_brejtau_interpolator->getNearestValue(input_values);
          mlog << MSG::DEBUG << " Rtau( "<<ptjet/1.e3<<" , "<<etajet<<" ) = "<< RejTau << endreq;
          if ( randnum < 1./RejTau ) 
            (*m_it)->setBTagged("SV1+IP3D");
        } else {
          double RejLight = m_interpolateBTagging ?
            m_brejpu_interpolator->interpolate(input_values) :
            m_brejpu_interpolator->getNearestValue(input_values);
          double dRb   = (*m_it)->getdRbquark();
          double dRc   = (*m_it)->getdRcquark();
          double dRtau = (*m_it)->getdRhadtau();
          if (dRb < 0.8 || dRc < 0.8 || dRtau < 0.8) {
            RejLight = m_interpolateBTagging ?
            m_brejnpu_interpolator->interpolate(input_values) :
            m_brejnpu_interpolator->getNearestValue(input_values);
          }
          mlog << MSG::DEBUG << " RLight( "<<ptjet/1.e3<<" , "<<etajet<<" ) = "<< RejLight << endreq;
          // correction factor only for light jet at zero order 
          // since the quality of reconstruction matters here 
          // (contrary to c jets, where physics dominates)
          if ( randnum < 1./(RejLight*m_correctionFactor) )
            (*m_it)->setBTagged("SV1+IP3D");
        }
      } //else{
        //(*m_it)->setLightTag();
      //}
    }
    return StatusCode::SUCCESS;
  }  
  
  StatusCode AtlfastB::atlfCal(){
    
    //It recalibrates the energy of all jets using parameterisation 
    //obtained with 'Z+jet' equivalent method
    
    MsgStream log( messageService(), name() ) ; 
    
    for(m_it=m_Jets->begin();m_it<m_Jets->end();++m_it){
      
      double ptjet=(*m_it)->pT();
      
      // Correction factor for b-jet momentum
      (*m_it)->setBTagCorrFactor("SV1+IP3D",fitcoreb(ptjet));
      
      // Correction factor for tau jets
      (*m_it)->setTauTagCorrFactor("Tau",1.);
      
      // Correction factor for b-jet momentum
      (*m_it)->setLightTagCorrFactor("Standard",fitcoreu(ptjet));      
      
    }
    
    return StatusCode::SUCCESS;

  }
  
  
  StatusCode AtlfastB::atlfTau(double epsitau){
    
    MsgStream log( messageService(), name() ) ; 
    //StatusCode sc;
    
    for(m_it=m_Jets->begin();m_it<m_Jets->end();++m_it){
      
      double randnum;
      
      double etajet=(*m_it)->eta();      
      double ptjet=(*m_it)->pT();
      int pdgid=(*m_it)->pdg_id();
      
      randnum=m_pRandFlatGenerator->shoot(m_pRandomEngine);
      
      deque<double> input_values;
      // Efficiencies are percentages in the input file
      input_values.push_back(epsitau*100.);
      // Efficiencies considered symmetric in eta
      input_values.push_back(fabs(etajet));
      // pT in GeV in input file
      input_values.push_back(ptjet/GeV);
      
      
      if(abs(pdgid)==15){
        
        // If the jet falls outside Interpolator limits, taueff = 0
        double taueff = m_taueff_interpolator->interpolate(input_values);
        
        if(randnum<taueff){
          
          (*m_it)->setTauTagged("Standard");
          log << MSG::DEBUG << "setTauTagged(\"Standard\"), pT = " << (*m_it)->pT() << endreq;
          
        }
        
      }else if(abs(etajet)<=2.5){
        
        double taurej = m_taurej_interpolator->interpolate(input_values);
        // taurej = 0 if values fall outside the Interpolator limits
        // Assume dummy rejection of 1e9
        if (!taurej) taurej = 1e9;
        
        if(randnum<(1./taurej)){
          (*m_it)->setTauTagged("Standard");
          log << MSG::DEBUG << "setTauTagged(\"Standard\"), pT = " << (*m_it)->pT() << endreq;
        }
        
      }
      
    }

    return StatusCode::SUCCESS;
    
  }
  
  
  StatusCode AtlfastB::tautag(double pt, 
                              double eta, 
                              double efftau, 
                              double &rjet, int &iflag){
    
    //It computes the corresponding jet rejection rjet
    //inputs: pt and eta of cluster, efftau-> eff(in %)  for tau identification
    
    MsgStream log( messageService(), name() ) ; 

    iflag=0;   


    double ptInGeV = pt/GeV;
    
    if((ptInGeV<15.)||(ptInGeV>150.)){
      
      
      log << MSG::DEBUG<<"Tau pt value out of range (15<pt<150 GeV): pt= "<<ptInGeV
          <<" GeV No tau tagging efficiencies will be applied"<<endreq;

      iflag=1;
      
    }
    if(abs(eta)>2.5){
      
      log << MSG::DEBUG<<"Tau eta value out of range (eta<2.5): eta= "<<eta<<" No tau tagging efficiencies will be applied"<<endreq;
      

      iflag=1;
      
    }
    

    //eta dependence
    double eff=efftau;
    double coefe1;
    double coefe3;
    double coef1;
    double coef2;
    
    if(abs(eta)<0.7){
      coefe1=1.35-0.0035*efftau;
      eff=efftau/coefe1;
    }
    if(abs(eta)>1.5){
      coefe3=0.70+0.0030*efftau;
      eff=efftau/coefe3;
    }
    if(eff>100){
      
      
      log << MSG::WARNING<<"Tau efficiency value > 100 ! eff= "<<eff<<endreq;
      iflag=1;
    }
    

    coef1=0.027+0.00024*ptInGeV;
    coef2=2.28+0.027*ptInGeV;

    coef1=0.027+0.00024*ptInGeV;
    coef2=2.28+0.027*ptInGeV;
    rjet=pow(10,(-coef1*eff+coef2));
    
  
    
    //     log << MSG::DEBUG<<"Tau rejection value computed in tautag= "<<rjet<<endreq;
    
    
    return StatusCode::SUCCESS;
    
  }
  

  StatusCode AtlfastB::atlfTau1P3P(){
    
    MsgStream log( messageService(), name() ) ; 
    //StatusCode sc;
    
    const double MaxPTinGeV = 150.0;
    
    for(m_it=m_Jets->begin();m_it<m_Jets->end();++m_it){
      double randnum;
      
      double etajet=(*m_it)->eta();      
      double ptjet=(*m_it)->pT();
      int pdgid=(*m_it)->pdg_id();
      
      randnum=m_pRandFlatGenerator->shoot(m_pRandomEngine);
      
      if(abs(pdgid)==15){

        if(ptjet/GeV > MaxPTinGeV){
        log << MSG::DEBUG<<"Jet with pT="<<ptjet/GeV<<" GeV; Outside range of Tau1P3P - Not Considered for tagging"<<endreq;
        }else{
        
          double tau1Peff = m_epsitau1P;
          double tau3Peff = m_epsitau3P;

          // Adjust tau3P efficiency due to lower efficiency at low PT
          if     ( ptjet/GeV >= 10.0 && ptjet/GeV < 20.0){ tau3Peff = 0.0; }
          else if( ptjet/GeV >= 20.0 && ptjet/GeV < 45.0){ tau3Peff = (0.0232*ptjet/GeV-0.044) * m_epsitau3P;}
          else if( ptjet/GeV >= 45.0 ){tau3Peff = m_epsitau3P;}

          if(randnum<tau1Peff){
            // tag as tau1P
            (*m_it)->setTauTagged("Tau1P3P:1prong");
            (*m_it)->setTauTagCorrFactor("Tau1P3P", 1. );
            log << MSG::DEBUG << "setTauTagged(\"Tau1P3P:1prong\"), pT = " << (*m_it)->pT() << endreq;
          }else if(randnum>=tau1Peff && randnum<(tau1Peff+tau3Peff) ){
            // tag as tau3P
            (*m_it)->setTauTagged("Tau1P3P:3prong");
            (*m_it)->setTauTagCorrFactor("Tau1P3P", 1. );
            log << MSG::DEBUG << "setTauTagged(\"Tau1P3P:3prong\"), pT = " << (*m_it)->pT() << endreq;
          }
        }       
      }else if(abs(etajet)<=2.5){

        if(ptjet/GeV > MaxPTinGeV){
        log << MSG::DEBUG<<"Jet with pT="<<ptjet/GeV<<" GeV; Outside range of Tau1P3P - Not Considered for tagging"<<endreq;
        }else{

          // Calculate probability of tagging jet as tau1P or tau3P
          double tau1Pprob = 0.0;
          double tau3Pprob = 0.0;

          if( m_epsitau1P != 0.0 ){
            double tau1Prej = -1.0;
            // Assign PT dependent rejection for Tau1P PT > 10 GeV
            if(ptjet/GeV >= 10. && ptjet/GeV < 20.){ 
              tau1Prej = exp(7.585-6.565*m_epsitau1P);}
            else if(ptjet/GeV >= 20. && ptjet/GeV < 30.){
              tau1Prej = exp(6.863-6.786*m_epsitau1P);}
            else if(ptjet/GeV >= 30. && ptjet/GeV < 40.){
              tau1Prej = exp(6.995-7.194*m_epsitau1P);}
            else if(ptjet/GeV >= 40. && ptjet/GeV < 50.){
              tau1Prej = exp(7.199-7.242*m_epsitau1P);}
            else if(ptjet/GeV >= 50. && ptjet/GeV < 60.){
              tau1Prej = exp(6.857-5.511*m_epsitau1P);}
            else if(ptjet/GeV >= 60.){
              tau1Prej = exp(7.344-6.331*m_epsitau1P);}

            if( tau1Prej != -1.0 ){ tau1Pprob = 1/tau1Prej; }
          }
        
          if( m_epsitau3P != 0.0 ){
            double tau3Prej = -1.0;
            // Assign PT dependent rejection for Tau3P PT > 20GeV
            if(ptjet/GeV >= 20. && ptjet/GeV < 30.){
              tau3Prej = exp(8.351-26.997*m_epsitau3P);}
            else if(ptjet/GeV >= 30. && ptjet/GeV < 40.){
              tau3Prej = exp(7.076-25.647*m_epsitau3P);}
            else if(ptjet/GeV >= 40. && ptjet/GeV < 50.){
              tau3Prej = exp(6.394-21.407*m_epsitau3P);}
            else if(ptjet/GeV >= 50. && ptjet/GeV < 60.){
              tau3Prej = exp(6.318-20.522*m_epsitau3P);}
            else if(ptjet/GeV >= 60.){
              tau3Prej = exp(6.477-19.238*m_epsitau3P);}

            if( tau3Prej != -1.0 ){ tau3Pprob = 1/tau3Prej; }
          } 
        
          if( randnum<tau1Pprob ){
            // Rescale Energy of fake Tau1P
            double corr = 0.0;
            double rand1P;
            // Do not let corr be small enough such that jet falls 
            // below 10 GeV jet reconstruction threshold.
            while( (corr * (*m_it)->pT())/GeV < 10.){
              // Get Random Number
              rand1P=m_pRandFlatGenerator->shoot(m_pRandomEngine);
              // Get PT bin of jet          
              int ptbin = 0;
              if(ptjet/GeV >= 10. && ptjet/GeV < 20.){      ptbin = 0; }
              else if(ptjet/GeV >= 20. && ptjet/GeV < 30.){ ptbin = 1; }
              else if(ptjet/GeV >= 30. && ptjet/GeV < 40.){ ptbin = 2; }
              else if(ptjet/GeV >= 40. && ptjet/GeV < 50.){ ptbin = 3; }
              else if(ptjet/GeV >= 50. && ptjet/GeV < 60.){ ptbin = 4; }
              else if(ptjet/GeV >= 60.){                    ptbin = 5; }
              // Get Scale Factor
              for(int ibin = 0; ibin<20; ibin++){
                if(m_corr_prob1P[ptbin][ibin] >= rand1P){
                  log << MSG::DEBUG << "pT = " << (*m_it)->pT() << " ptbin = " << ptbin << " rand1P = " << rand1P << " corr = " << m_corr1P3P[ibin] << endreq;
                  corr = m_corr1P3P[ibin];
                  break;
                }
              } 
            }
            // Tag as Tau1P and rescale momentum.
            log << MSG::DEBUG<<"Jet tagged as Tau1P. Initial Jet PT = "<<ptjet<<"; After rescaling PT = "<<corr*(*m_it)->pT()/GeV <<endreq;
            (*m_it)->setTauTagged("Tau1P3P:1prong");
            (*m_it)->setTauTagCorrFactor("Tau1P3P",corr);
                  
          }else if( randnum>=tau1Pprob && randnum<(tau1Pprob+tau3Pprob) ){
          
            // Rescale Energy of fake Tau3P
            double corr = 0.0;
            double rand3P;
            // Do not let corr be small enough such that jet falls 
            // below 10 GeV jet reconstruction threshold.
            while( (corr * (*m_it)->pT())/GeV < 10.){
              // Get Random Number
              rand3P=m_pRandFlatGenerator->shoot(m_pRandomEngine);
              // Get PT bin of jet          
              int ptbin = 0;
              if(ptjet/GeV >= 10. && ptjet/GeV < 20.){      ptbin = 0; }
              else if(ptjet/GeV >= 20. && ptjet/GeV < 30.){ ptbin = 1; }
              else if(ptjet/GeV >= 30. && ptjet/GeV < 40.){ ptbin = 2; }
              else if(ptjet/GeV >= 40. && ptjet/GeV < 50.){ ptbin = 3; }
              else if(ptjet/GeV >= 50. && ptjet/GeV < 60.){ ptbin = 4; }
              else if(ptjet/GeV >= 60.){                    ptbin = 5; }
              // Get Scale Factor
              for(int ibin = 0; ibin<20; ibin++){
                if(m_corr_prob3P[ptbin][ibin] >= rand3P){
                  log << MSG::DEBUG << "pT = " << (*m_it)->pT() << " ptbin = " << ptbin << " rand3P = " << rand3P << " corr = " << m_corr1P3P[ibin] << endreq;
                  corr = m_corr1P3P[ibin];
                  break;
                }
              }
            }
            if( corr * (*m_it)->pT()/GeV >= 20.0 ){
              // If rescaled ET >= 20 GeV, tag as Tau3P and rescale momentum.
              // If rescaled ET <  20 GeV, do not tag (tau3P low limit is 20GeV)
              log << MSG::DEBUG<<"Jet tagged as Tau3P. Initial Jet PT = "<<ptjet<<"; After rescaling PT = "<<corr*(*m_it)->pT()/GeV <<endreq;
              (*m_it)->setTauTagged("Tau1P3P:3prong");
              (*m_it)->setTauTagCorrFactor("Tau1P3P",corr);
            }
          }
        }
      }
    }
    return StatusCode::SUCCESS;
  }

  
  StatusCode  AtlfastB::atlfTauVeto(int ind){
    
    MsgStream log( messageService(), name() ) ; 
    
    
    StatusCode sc;
    
    for(m_it=m_Jets->begin();m_it<m_Jets->end();++m_it){
      
      double randnum;    
      double efftau;
      double effjet;
      
      double ptjet=(*m_it)->pT();
      int pdgid=(*m_it)->pdg_id();

      sc=tauveto(ptjet,ind,efftau,effjet);
      if(sc.isFailure()){
        log << MSG::ERROR<<"Cannot compute tau jet rejection"<<endreq;
        return StatusCode::FAILURE;
      }
      
      randnum=m_pRandFlatGenerator->shoot(m_pRandomEngine);
      

      if(abs(pdgid)==15) {
        if(randnum<efftau/100.){
          // Need to unset the tau tag here
          //(*m_it)->setLightTag();
        }
      }else{
        if(randnum>effjet/100.){
          (*m_it)->setTauTagged("Veto");
        }
      }
    }
    
    return StatusCode::SUCCESS;
  }
  
  
  
  StatusCode AtlfastB::tauveto(double pt, int ind, double &efftau, double &effjet){ 
    MsgStream log( messageService(), name() ) ; 
    
    

    //input: pt of cluster
    //ind=1->fix efftau to 5% and gives effjet
    //ind=2->fix effjet to 90% and gives efftau
    //output: efftau tau efficiency in % or effjet jets efficiency in %
    
    
    //using calo criteria + tracks
    
    
    double ptInGeV = pt/GeV;

    if(ind==1){
     
      efftau=5.;
      effjet=36.+1.5*ptInGeV-0.01*pow(ptInGeV,2);
      if(ptInGeV>60.) effjet=90.;
    }
    else if(ind==2){
      effjet=90.;
      efftau=27.-0.35*ptInGeV;
      if(ptInGeV>60.) efftau=5.;

    }else{
      
      log << MSG::ERROR<<"Wrong ind values for tauveto method. ind= "<<ind<<"ind=1 or ind=2"<<endreq;
      
      return StatusCode::FAILURE;
    }
    
    
    
    return StatusCode::SUCCESS;
  }
  
  StatusCode AtlfastB::atlfTrigMuo(){
    MsgStream log( messageService(), name() ) ; 
    
    log << MSG::INFO<<"atlfTrigMuo is not yet implemented!!!!!!"<<endreq;
    
    
    return StatusCode::SUCCESS;
  }
  double AtlfastB::fitcoreb(double pt){
    
    
    MsgStream log( messageService(), name() ) ; 
    
    double core;
    double a0,a1,a2,a3,a4,a5;
    double xc;
    
    
    double ptInGeV = pt/GeV;
    
    if(ptInGeV<10.){
      
      core=0.;
    }else if(ptInGeV<55.){
      
      a0=1.26694; // Changed from 1.2715 to fix discontinuity
      a1=0.12241;
      a2=-0.10480e-01;
      a3=0.33310e-03;
      a4=-0.47454e-05;
      a5=0.25436e-07;
      core=a0+a1*ptInGeV+a2*pow(ptInGeV,2)+
        a3*pow(ptInGeV,3)+a4*pow(ptInGeV,4)+
        a5*pow(ptInGeV,5);
      core=core*1.006;
    }else if(ptInGeV<200.){
      
      a0=1.18;
      a1=-0.16672e-02;
      a2=0.44414e-05;
      core=a0+a1*ptInGeV+a2*pow(ptInGeV,2);
      
    }else{
      
      xc=200.;
      a0=1.18;
      a1=-0.16672e-02;
      a2=0.44414e-05;
      core=a0+a1*xc+a2*pow(xc,2);
    }
    
    // core corrects pt**1
    
    return core;
    
  }
  
  double AtlfastB::fitcoreu(double pt){
    
    MsgStream log( messageService(), name() ) ;     
    
    double core;
    double a0,a1,a2,a3,a4,a5;
    double xc;
    
    double ptInGeV = pt/GeV;
    
    if(ptInGeV<10.){
      
      core=0.;
      
      
      
      
      
    }else if((ptInGeV<45.)&&(ptInGeV>10.)){
      
      a0=1.5065; // Changed from 1.5085 to fix discontinuity
      a1=0.31468e-01;
      a2=-0.36973e-02;
      a3=0.11220e-03;
      a4=-0.13921e-05;
      a5=0.61538e-08;
      core=a0+a1*ptInGeV+a2*pow(ptInGeV,2)+
        a3*pow(ptInGeV,3)+a4*pow(ptInGeV,4)+
        a5*pow(ptInGeV,5);
      
    }else if(ptInGeV<200.){
      
      a0=1.18;
      a1=-0.16672e-02;
      a2=0.44414e-05;
      core=a0+a1*ptInGeV+a2*pow(ptInGeV,2);
      core=core/1.025;
    }else{
      
      xc=200.;
      a0=1.18;
      a1=-0.16672e-02;
      a2=0.44414e-05;
      core=a0+a1*xc+a2*pow(xc,2);
      core=core/1.025;
    }
    
    
    // core corrects pt**1
    return core;
    
  }

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

---