// ================================================
 // 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"
 
 
 /** @brief Applies jet tagging efficiencies and momentum correction factors.
    *
    * The original jet flavour labels from JetMaker are rearranged here to
    * simulate the overall efficiencies measured in full simulation.
    * Correction factors are also computed to match different tagging hypotheses.
    *
    * @image html AtlfastB.jpg "AtlfastB sequence"
    */
 
 //--------------------------------
 // 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;
     /** Sets the b-tagging regime to be used. */
     m_atlfBNSet           =  10;
     //!< tau-jet rejection. @see m_atlfBNSet
     m_epsitau             =  0.5;
     //!< Flag to configure tauveto. @see tauveto(double pt, int ind, double &efftau, double &effjet)
     m_indtauveto          =  1;
     //!< input file for light jet correction factors
     m_corrjfile           =  "atlfastDatafiles/AtlfastBjet.dat";
     //!< input file for c-jet correction factors
     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) { // use TDR parameterisation?
       //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
 

Due to the LXR bug, the updates fail sometimes to remove references to deleted files. The Saturday's full rebuilds fix these problems
This page was automatically generated by the LXR engine.