NtupleReader/NtupleReader/MyTauJet.ixx

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// 24.10.2006, AUTHOR: THIES EHRICH
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
//:: IMPLEMENTATION OF INLINE METHODS DEFINED IN THE CLASS MyTauJet ::
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

//*****************************************************************************

  //:::::::::::::::::
  //:: METHOD init ::
  //:::::::::::::::::

inline void MyTauJet::init(const TLorentzVector & p,
                           const double & charge, 
                           const int & rec_flag,
                           const int & index,
                           const double & emRadius,
                           const double & isolationFraction,
                           const double & centralityFraction,
                           const double & stripWidth2,
                           const double & nStripCells,
                           const double & llh,
                           const double & lowPtTauEleDiscriminant,
                           const double & tauENeuralNetwork,
                           const double & tauJetNeuralnetwork,
                           const double & etHadCalib, 
                           const double & etEMCalib,
                           const int & nb_tracks,
                           const TLorentzVector & track_1, 
                           const TLorentzVector & track_2, 
                           const TLorentzVector & track_3,
                           const int & track_TRTHits_1,
                           const int & track_TRTHighThresholdHits_1,
                           const int & track_TRTHits_2,
                           const int & track_TRTHighThresholdHits_2,
                           const int & track_TRTHits_3,
                           const int & track_TRTHighThresholdHits_3
                           ) {

    m_dataset_info = &MyDatasetInfo::Instance();
    m_tlv = p;
    m_charge = charge;
    m_rec_flag = rec_flag;
    m_index = index;
    m_emRadius= emRadius;
    m_isolationFraction=isolationFraction;
    m_centralityFraction=centralityFraction;
    m_stripWidth2=stripWidth2;
    m_nStripCells=nStripCells;
    m_llh = llh;
    m_lowPtTauEleDiscriminant = lowPtTauEleDiscriminant ;
    m_tauENeuralNetwork = tauENeuralNetwork;
    m_tauJetNeuralnetwork = tauJetNeuralnetwork;
    m_etHadCalib = etHadCalib;
    m_etEMCalib = etEMCalib;
    m_nb_tracks = nb_tracks;
    m_track_1 = track_1;
    m_track_2 = track_2;
    m_track_3 = track_3;
    m_track_TRTHits_1 = track_TRTHits_1;
    m_track_TRTHighThresholdHits_1 = track_TRTHighThresholdHits_1;
    m_track_TRTHits_2 = track_TRTHits_2;
    m_track_TRTHighThresholdHits_2 = track_TRTHighThresholdHits_2;
    m_track_TRTHits_3 = track_TRTHits_3;
    m_track_TRTHighThresholdHits_3 = track_TRTHighThresholdHits_3;
    m_type = "tau_jet";
    m_pdgId = TMath::Sign(1,(int)m_charge)*-15;
    m_matchingflag = 0;

    return;

  }

//*****************************************************************************

  //::::::::::::::::::::::::::::::::
  //:: METHOD reconstruction_flag ::
  //::::::::::::::::::::::::::::::::

inline int MyTauJet::reconstruction_flag(void) const {

    return m_rec_flag;

  }

//*****************************************************************************

  //::::::::::::::::
  //:: METHOD llh ::
  //::::::::::::::::

inline double MyTauJet::llh(void) const {

    return m_llh;

  }

//*****************************************************************************

  //::::::::::::::::::::::::::::::::::::
  //:: METHOD lowPtTauEleDiscriminant ::
  //::::::::::::::::::::::::::::::::::::

inline double MyTauJet::lowPtTauEleDiscriminant(void) const {

    return m_lowPtTauEleDiscriminant;

  }
//*****************************************************************************

  //::::::::::::::::::::::::::::::
  //:: METHOD tauENeuralNetwork ::
  //::::::::::::::::::::::::::::::

inline double MyTauJet::tauENeuralNetwork(void) const {

    return m_tauENeuralNetwork;

  }
//*****************************************************************************

  //::::::::::::::::::::::::::::::::
  //:: METHOD tauJetNeuralnetwork ::
  //::::::::::::::::::::::::::::::::

inline double MyTauJet::tauJetNeuralnetwork(void) const {

    return m_tauJetNeuralnetwork;

  }

//*****************************************************************************

  //::::::::::::::::::::::::::::::::::::
  //:: METHOD set_reconstruction_flag ::
  //::::::::::::::::::::::::::::::::::::

inline void MyTauJet::set_reconstruction_flag(const int & flag) {

    m_rec_flag = flag;
    return;

  }

//*****************************************************************************

  //:::::::::::::::::::::::
  //:: METHOD emRadius ::
  //:::::::::::::::::::::::

inline double MyTauJet::emRadius(void) const{

    return m_emRadius;

  }
//*****************************************************************************

  //:::::::::::::::::::::::
  //:: METHOD isolationFraction ::
  //:::::::::::::::::::::::

inline double MyTauJet::isolationFraction(void) const{

    return m_isolationFraction;

  }
//*****************************************************************************

  //:::::::::::::::::::::::
  //:: METHOD etHadCalib ::
  //:::::::::::::::::::::::

inline double MyTauJet::centralityFraction(void) const{

    return m_centralityFraction;

  }
//*****************************************************************************

  //:::::::::::::::::::::::
  //:: METHOD etHadCalib ::
  //:::::::::::::::::::::::

inline double MyTauJet::stripWidth2(void) const{

    return m_stripWidth2;

  }
//*****************************************************************************

  //:::::::::::::::::::::::
  //:: METHOD nStripCells ::
  //:::::::::::::::::::::::

inline double MyTauJet::nStripCells(void) const{

    return m_nStripCells;

  }
//*****************************************************************************

  //:::::::::::::::::::::::
  //:: METHOD etHadCalib ::
  //:::::::::::::::::::::::

inline double MyTauJet::ethadcalib(void) const{

    return m_etHadCalib;

  }
//*****************************************************************************

  //::::::::::::::::::::::
  //:: METHOD etEMCalib ::
  //::::::::::::::::::::::

inline double MyTauJet::etemcalib(void) const {

    return m_etEMCalib;

  }
//*****************************************************************************

  //::::::::::::::::::::::
  //:: METHOD nb_tracks ::
  //::::::::::::::::::::::

inline int MyTauJet::nb_tracks(void) const {

    return m_nb_tracks;

  }
//*****************************************************************************

  //:::::::::::::::::::::::
  //:: METHOD pt_track_1 ::
  //:::::::::::::::::::::::

inline TLorentzVector MyTauJet::track_1(void) const {

    return m_track_1;

  }
//*****************************************************************************

  //:::::::::::::::::::::::
  //:: METHOD pt_track_2 ::
  //:::::::::::::::::::::::

inline TLorentzVector MyTauJet::track_2(void) const {

    return m_track_2;

  }
//*****************************************************************************

  //:::::::::::::::::::::::
  //:: METHOD pt_track_3 ::
  //:::::::::::::::::::::::

inline TLorentzVector MyTauJet::track_3(void) const {

    return m_track_3;

  }
//*****************************************************************************

  //::::::::::::::::::::::::::::
  //:: METHOD track_TRTHits_1 ::
  //::::::::::::::::::::::::::::

inline int MyTauJet::track_TRTHits_1(void) const {

    return m_track_TRTHits_1;

  }
//*****************************************************************************

  //:::::::::::::::::::::::::::::::::::::::::
  //:: METHOD track_TRTHighThresholdHits_1 ::
  //:::::::::::::::::::::::::::::::::::::::::

inline int MyTauJet::track_TRTHighThresholdHits_1(void) const {

    return m_track_TRTHighThresholdHits_1;

  }
//*****************************************************************************

  //::::::::::::::::::::::::::::
  //:: METHOD track_TRTHits_2 ::
  //::::::::::::::::::::::::::::

inline int MyTauJet::track_TRTHits_2(void) const {

    return m_track_TRTHits_2;

  }
//*****************************************************************************

  //:::::::::::::::::::::::::::::::::::::::::
  //:: METHOD track_TRTHighThresholdHits_2 ::
  //:::::::::::::::::::::::::::::::::::::::::

inline int MyTauJet::track_TRTHighThresholdHits_2(void) const {

    return m_track_TRTHighThresholdHits_2;

  }
//*****************************************************************************

  //::::::::::::::::::::::::::::
  //:: METHOD track_TRTHits_3 ::
  //::::::::::::::::::::::::::::

inline int MyTauJet::track_TRTHits_3(void) const {

    return m_track_TRTHits_3;

  }
//*****************************************************************************

  //:::::::::::::::::::::::::::::::::::::::::
  //:: METHOD track_TRTHighThresholdHits_3 ::
  //:::::::::::::::::::::::::::::::::::::::::

inline int MyTauJet::track_TRTHighThresholdHits_3(void) const {

    return m_track_TRTHighThresholdHits_3;

  }
//*****************************************************************************

  //::::::::::::::::::
  //:: METHOD clone ::
  //::::::::::::::::::

inline MyTauJet* MyTauJet::Clone(void) const {

    double Px = this->tlv().Px();
    double Py = this->tlv().Py();
    double Pz = this->tlv().Pz();
    double E = this->tlv().E();
    double charge = this->charge();
    int rec_flag = this->reconstruction_flag();
    int index = this->index();
    double llh = this->llh();
    double emRadius = this->emRadius();
    double centralityFraction= this->centralityFraction();
    double isolationFraction= this->isolationFraction();
    double stripWidth2= this->stripWidth2();
    double nStripCells= this->nStripCells();
    double lowPtTauEleDiscriminant = this->lowPtTauEleDiscriminant();
    double tauENeuralNetwork = this->tauENeuralNetwork();
    double tauJetNeuralnetwork = this->tauJetNeuralnetwork();
    double etHadCalib = this->ethadcalib();
    double etEMCalib = this->etemcalib();
    int nb_tracks = this->nb_tracks();
    double track_px_1 = this->track_1().Px();
    double track_py_1 = this->track_1().Py();
    double track_pz_1 = this->track_1().Pz();
    double track_e_1 = this->track_1().E();
    double track_px_2 = this->track_2().Px();
    double track_py_2 = this->track_2().Py();
    double track_pz_2 = this->track_2().Pz();
    double track_e_2 = this->track_2().E();
    double track_px_3 = this->track_2().Px();
    double track_py_3 = this->track_3().Py();
    double track_pz_3 = this->track_3().Pz();
    double track_e_3 = this->track_3().E();
    int track_TRTHits_1 = this->track_TRTHits_1();
    int track_TRTHighThresholdHits_1 = this->track_TRTHighThresholdHits_1();
    int track_TRTHits_2 = this->track_TRTHits_2();
    int track_TRTHighThresholdHits_2 = this->track_TRTHighThresholdHits_2();
    int track_TRTHits_3 = this->track_TRTHits_3();
    int track_TRTHighThresholdHits_3 = this->track_TRTHighThresholdHits_3();
    


    MyTauJet* cloned_taujet 
        = new MyTauJet(TLorentzVector(Px,Py,Pz,E),charge,
                       rec_flag, index,
                       emRadius,
                       centralityFraction,
                       isolationFraction,
                       stripWidth2,
                       nStripCells,
                       llh, lowPtTauEleDiscriminant, 
                       tauENeuralNetwork, tauJetNeuralnetwork,
                       etHadCalib, etEMCalib, nb_tracks,
                       TLorentzVector(track_px_1,track_py_1,track_pz_1, track_e_1), 
                       TLorentzVector(track_px_2,track_py_2,track_pz_2, track_e_2),
                       TLorentzVector(track_px_3,track_py_3,track_pz_3, track_e_3),
                       track_TRTHits_1,
                       track_TRTHighThresholdHits_1,
                       track_TRTHits_2,
                       track_TRTHighThresholdHits_2,
                       track_TRTHits_3,
                       track_TRTHighThresholdHits_3
                       ); 
    
    return cloned_taujet;

  }
//*****************************************************************************

  //:::::::::::::::::::::::::::::
  //:: METHOD get corrected Et ::
  //:::::::::::::::::::::::::::::

inline TLorentzVector MyTauJet::tlv_corr(void) const {

    double m_gev;
    m_gev = m_dataset_info->get_scale();
    
    double m_ptetacorrectionsntr1[9] = { 0.8339,   0.8125,   0.8254,
                                         0.8747,   0.8583,   0.8594,
                                         0.9088,   0.9013,   0.8922};
    double m_ptetacorrectionsntr23[9]= { 0.9000,   0.8593,   0.9034,
                                         0.9208,   0.8791,   0.9132,
                                         0.9359,   0.9231,   0.9033};

    const int m_nptbins  = 3;
    const int m_netabins = 3;

    double m_ptpoints[m_nptbins+1] = {15000.*m_gev, 35000.*m_gev, 150000.*m_gev, 0.};
    double m_etapoints[m_netabins+1] = {0.25, 1.0, 2.0, 0.};

    double et_corr = m_etEMCalib+m_etHadCalib;
    double eta = fabs( m_tlv.Eta() );
    double corr = 1, fudge = 1.011;

    unsigned int lowpt_idx = 0;
    unsigned int loweta_idx = 0;

    double lowpt_frac = 0;
    double loweta_frac = 0;
                        
    double e_corr;
    double p_corr;
    double pt_corr;
    TLorentzVector m_tlv_corr;
                        
    while ( lowpt_idx < m_nptbins-1 && et_corr > m_ptpoints[lowpt_idx+1] )
        lowpt_idx++;

    lowpt_frac = ( m_ptpoints[lowpt_idx+1] - et_corr ) / ( m_ptpoints[lowpt_idx+1] -
                                                           m_ptpoints[lowpt_idx] );
    // will be >1 if et_corr is out of bounds

    if ( lowpt_frac > 1 )
        lowpt_frac =  1;
    if( lowpt_frac < 0 ) {   
        //should never happen, only if ptNumberOfBins is set wrong (which is checked now)
        std::cerr << "FIXME: lowpt_frac < 0 !!" << std::endl;
    }

    while ( loweta_idx < m_netabins-1 && eta > m_etapoints[loweta_idx+1] )
        loweta_idx++;

    loweta_frac = ( m_etapoints[loweta_idx+1] - eta ) / 
        ( m_etapoints[loweta_idx+1] - m_etapoints[loweta_idx] );
    // will be >1 if eta is out of bounds, 

    if ( loweta_frac > 1)
        loweta_frac = 1;
    if( loweta_frac < 0 ) {    //should never happen, only if etaNumberOfBins is set wrong
        std::cerr << "FIXME: loweta_frac < 0 !!" << std::endl;
    }

    double coeff_matrix[2][2] = { {0, 0}, {0, 0} };

    if ( m_nb_tracks <= 1 ) {
        coeff_matrix[0][0] = m_ptetacorrectionsntr1[lowpt_idx*m_netabins+loweta_idx];
        if( lowpt_idx < m_nptbins-1 )
            coeff_matrix[1][0] = m_ptetacorrectionsntr1[(lowpt_idx+1)*m_netabins+loweta_idx];
        if( loweta_idx < m_netabins-1 )
            coeff_matrix[0][1] = m_ptetacorrectionsntr1[lowpt_idx*m_netabins+(loweta_idx+1)];
        if( lowpt_idx < m_nptbins-1 && loweta_idx < m_netabins-1 )
            coeff_matrix[1][1] = m_ptetacorrectionsntr1[(lowpt_idx+1)*m_netabins+(loweta_idx+1)];
    } else {
        coeff_matrix[0][0] = m_ptetacorrectionsntr23[lowpt_idx*m_netabins+loweta_idx];
        if( lowpt_idx < m_nptbins-1 )
            coeff_matrix[1][0] = m_ptetacorrectionsntr23[(lowpt_idx+1)*m_netabins+loweta_idx];
        if( loweta_idx < m_netabins-1 )
            coeff_matrix[0][1] = m_ptetacorrectionsntr23[lowpt_idx*m_netabins+(loweta_idx+1)];
        if( lowpt_idx < m_nptbins-1 && loweta_idx < m_netabins-1 )
            coeff_matrix[1][1] = m_ptetacorrectionsntr23[(lowpt_idx+1)*m_netabins+(loweta_idx+1)];
    }

    corr = ( coeff_matrix[0][0] * lowpt_frac * loweta_frac ) 
        + ( coeff_matrix[1][0] * (1-lowpt_frac) * loweta_frac );
    
    corr += ( coeff_matrix[0][1] * lowpt_frac * (1-loweta_frac) ) 
        + ( coeff_matrix[1][1] * (1-lowpt_frac) * (1-loweta_frac) );

    et_corr *= corr*fudge;
                        
    //create lorentzvector//
    
    e_corr=sqrt(pow(et_corr,2)+pow(et_corr,2)*pow(cos(m_tlv.Theta())/sin(m_tlv.Theta()),2));
    
    p_corr=sqrt(pow(et_corr,2)+pow(et_corr,2)*pow(cos(m_tlv.Theta())/sin(m_tlv.Theta()),2)
                -pow(1776.99*m_gev,2));
    
    pt_corr=p_corr*sin(m_tlv.Theta());
    
    m_tlv_corr. SetPtEtaPhiE(pt_corr,m_tlv.Eta(),m_tlv.Phi(),e_corr);

                        
    return m_tlv_corr;
  }

//*****************************************************************************

  //::::::::::::::::::::::::::
  //:: METHOD PrintParticle ::
  //::::::::::::::::::::::::::

inline void MyTauJet::PrintParticle(std::string option) {

    std::cout << "   rec_flag, emRadius, isolFrac: " 
              << Form("%i, %.2f, %.2f",
                      this->reconstruction_flag(),
                      this->emRadius(),
                      this->isolationFraction() ) << std::endl;
              
    if(option=="v"){
        
        std::cout << "   centrFrac, strWi2, nStrCells: " 
                  << Form("%.2f, %.2f, %.2f",
                          this->centralityFraction(),
                          this->stripWidth2(),
                          this->nStripCells() ) << std::endl;

        std::cout << "   lowPtTauEleDiscr, tauENN:     " 
                  << Form("%.2f, %.2f",
                          m_lowPtTauEleDiscriminant,
                          m_tauENeuralNetwork) << std::endl;

        std::cout << "   tauJetNeuralnetwork, llh:     " 
                  << Form("%.2f, %.2f",
                          m_tauJetNeuralnetwork,
                          m_llh ) << std::endl;

        std::cout << "   ethadcalib, etemcalib, n_trk: " 
                  << Form("%.2f, %.2f, %i",
                          m_etHadCalib,
                          m_etEMCalib,
                          m_nb_tracks ) << std::endl;
    
        if(this->nb_tracks()>0){
            std::cout << "   trk1-TRT HiThrHits, Hits:     " 
                      << Form("%i, %i",
                              m_track_TRTHighThresholdHits_1,
                              m_track_TRTHits_1 ) << std::endl;
            std::cout << "   trk1 Pt, Eta, Phi:            "
                      << Form("%.2f, %.2f, %.2f",
                              m_track_1.Pt(), 
                              m_track_1.Eta(), 
                              m_track_1.Phi()) << std::endl;
        }
    
        if(this->nb_tracks()>1){
            std::cout << "   trk2-TRT HiThrHits, Hits:     " 
                      << Form("%i, %i",
                              m_track_TRTHighThresholdHits_2,
                              m_track_TRTHits_2 ) << std::endl;
            std::cout << "   trk2 Pt, Eta, Phi:            "
                      << Form("%.2f, %.2f, %.2f",
                              m_track_2.Pt(), 
                              m_track_2.Eta(), 
                              m_track_2.Phi()) << std::endl;
        }
    
        if(this->nb_tracks()>2){
            std::cout << "   trk2-TRT HiThrHits, Hits:     " 
                      << Form("%i, %i",
                              m_track_TRTHighThresholdHits_3,
                              m_track_TRTHits_3 ) << std::endl;
            std::cout << "   trk3 Pt, Eta, Phi:            "
                      << Form("%.2f, %.2f, %.2f",
                              m_track_3.Pt(), 
                              m_track_3.Eta(), 
                              m_track_3.Phi()) << std::endl;
        }

    }

    return;
  }

---