NtupleReader/src/MyTools.cxx

//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// 10.12.2006, AUTHOR: ALL
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
//:: IMPLEMENTATION OF METHODS DEFINED IN THE CLASS MyTools ::
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::

//::::::::::::::::::
//:: HEADER FILES ::
//::::::::::::::::::

#include "NtupleReader/MyTools.h"


//::::::::::::::::::::::::
//:: NAMESPACE SETTINGS ::
//::::::::::::::::::::::::

using namespace std;


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

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

void MyTools::init() {

    m_nb_multiple_matches = 0;
}

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

//:::::::::::::::::::::
//:: METHOD destruct ::
//:::::::::::::::::::::

void MyTools::destruct(void) {

    return;

}



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

//:::::::::::::::::::::::::::::::::::
//:: METHOD find_matching_particle ::
//:::::::::::::::::::::::::::::::::::

MyParticle*
MyTools::find_matching_particle( TLorentzVector particle_tlv, 
                                 std::vector<MyParticle*> & particle_vector, 
                                 double deltaR_cut,
                                 bool bool_matchingflag,
                                 bool bool_usematching){

    // Variables //

    double deltaR = 0.;
    double deltaR_min = 10.5;

    int k_min = -1;
    int nb_foundparticles =  0;
    //double deltaEta = 0;
    //double deltaPhi = 0;

    // loop over particles in vector // 

    for (unsigned int k=0; k<particle_vector.size(); k++) {

        // don't use particle, if it has already been matched before // 

        if(bool_usematching) {
            if( !(particle_vector[k]->matchingflag() == 0) )  continue; 
        }

        //deltaEta = particle->tlv().Eta()-particle_vector[k]->tlv().Eta();
        //deltaPhi = particle->tlv().Phi()-particle_vector[k]->tlv().Phi();
        //if(deltaPhi > 3.14159) deltaPhi = deltaPhi- 2.*3.1415; 
        //if(deltaPhi < -3.14159) deltaPhi = deltaPhi+ 2.*3.1415; 
        //deltaR = sqrt(pow(deltaEta,2)+pow(deltaPhi,2));

        deltaR = particle_tlv.DeltaR( particle_vector[k]->tlv() );
        //check the influence of a delta pt cut in matching
        //double delta_pt_rel;
        //delta_pt_rel = TMath::Abs((particle->tlv().Pt() - particle_vector[k]->tlv().Pt())) / particle->tlv().Pt();
        //if(delta_pt_rel<0.1){

        // select particle with dR < dR_cut                               //
        // if more than one canditate select the one with the smallest dR //

        if( deltaR < deltaR_cut && deltaR < deltaR_min){
            k_min = k;
            deltaR_min = deltaR; 
            nb_foundparticles++;
        }
        //} // end of influence of a delta pt cut

    }

    // if more than one particle with dR < dR_cut was found //
    // increase m_nb_multiple_matches                       //

    if(nb_foundparticles>1){ 
        m_nb_multiple_matches++;
    }

    if(nb_foundparticles==0){
        return NULL;
    }
    else{
        if(bool_matchingflag) particle_vector[k_min]->set_matchingflag(1);
        return particle_vector[k_min];
    }

  }

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

//:::::::::::::::::::::::::::::::::::
//:: METHOD find_matching_particle ::
//:::::::::::::::::::::::::::::::::::

MyParticle*
MyTools::find_matching_particle( MyParticle* particle, 
                                 std::vector<MyParticle*> & particle_vector, 
                                 double deltaR_cut,
                                 bool bool_matchingflag,
                                 bool bool_usematching){
    return find_matching_particle(particle->tlv(),
                                 particle_vector, 
                                 deltaR_cut,
                                 bool_matchingflag,
                                 bool_usematching);
}

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

Int_t 
MyTools::nb_multiple_matches(void){

    return m_nb_multiple_matches;

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

void 
MyTools::clear_vector(std::vector<MyParticle*> & particle_vector){

    for(unsigned int k =0;k<particle_vector.size(); k++ ){
        particle_vector[k]->set_matchingflag(0);
    }

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

void 
MyTools::clear_nb_multiple_matches(){

    m_nb_multiple_matches = 0;

}

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

//:::::::::::::::::::::::::::
//:: METHOD remove_overlap ::
//:::::::::::::::::::::::::::

std::vector<MyParticle*>
MyTools::remove_overlap( std::vector<MyParticle*> vec_part_tb_rm,
                         std::vector<MyParticle*> vec_part_rm_with, 
                         double deltaR_cut) {

    std::vector<MyParticle*> vec_tmp;

    MyParticle* tmp_part_overlaps = NULL;

    // loop over given vector
    for (unsigned int k=0; k<vec_part_tb_rm.size(); k++){

        tmp_part_overlaps = 
            find_matching_particle( vec_part_tb_rm[k], vec_part_rm_with, 
                                    deltaR_cut, false, false);

        // if maching particle was found go to the next in the vector
        if(tmp_part_overlaps!=NULL){
            continue;
        }

        vec_tmp.push_back( vec_part_tb_rm[k] );
    }

    return vec_tmp;

}

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

//:::::::::::::::::::::::::::
//:: METHOD sort_particles ::
//:::::::::::::::::::::::::::

std::vector<MyParticle*>
MyTools::sort_particles( std::vector<MyParticle*> p_vector ){


    std::vector<MyParticle*> sorted_p_vector(0);

    unsigned int orig_size = p_vector.size();

    for (unsigned int cycle=0; cycle<orig_size; cycle++){

        int    k_max=-1; 
        double abspt_max=-1;

        for (unsigned int k=0; k<p_vector.size(); k++) {

            double abspt =  TMath::Abs(p_vector[k]->tlv().Pt());

            if( abspt > abspt_max ){
                k_max = k;
                abspt_max = abspt; 
            }

        }

        //cout << "in sorter, cycle: " << cycle 
        //   << ", abspt_max: " << abspt_max << endl; 

        sorted_p_vector.push_back( p_vector[k_max] );
        p_vector.erase(p_vector.begin()+k_max);

    }

    return sorted_p_vector;

}

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

//::::::::::::::::::::::::
//:: METHOD check_order ::
//::::::::::::::::::::::::

void
MyTools::check_order( std::vector<MyParticle*> p_vector ){

    double abspt_min=0;

    if (p_vector.size()==0) return;

    for (unsigned int k=0; k<p_vector.size(); k++) {

        double abspt =  TMath::Abs(p_vector[k]->tlv().Pt());

        if(k!=0 && abspt>abspt_min){
            //cout << "k: " << k << "abspt_min: " << abspt_min << endl; 
            //cout << "k: " << k << "abspt    : " << abspt  << endl; 
            cerr << endl
                 << "Class MyTools, method check_order: ERROR!\n"
                 << "vector: " << p_vector[0]->type() << " not sorted" << endl;
            //exit(1);
        }

        //cout << "k: " << k << "abspt_min: " << abspt_min << endl; 
        abspt_min = abspt; 
    }

    return; 
}


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

//:::::::::::::::::::::::
//:: METHOD pt_eta_cut ::
//:::::::::::::::::::::::

std::vector<MyParticle*> 
MyTools::pt_eta_cut( const vector<MyParticle*> p_vector,
                     const double pt_min,  const double pt_max,
                     const double eta_min, const double eta_max) {

    std::vector<MyParticle*> tmp_v;
    tmp_v.clear();

    for (unsigned int i=0; i<p_vector.size(); i++){

        double pt  = (p_vector[i])->tlv().Pt();
        double eta = TMath::Abs((p_vector[i])->tlv().Eta());

        if(pt_min>-1){
            if(pt < pt_min) continue; 
        }
        if(pt_max>-1){
            if(pt > pt_max) continue; 
        }
        if(eta_min>-1){
            if(eta < eta_min) continue; 
        }
        if(eta_max>-1){
            if(eta > eta_max) continue; 
        }

        tmp_v.push_back(p_vector[i]);

    }

    return tmp_v;
}

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

//:::::::::::::::::::::::
//:: METHOD et_eta_cut ::
//:::::::::::::::::::::::

std::vector<MyParticle*> 
MyTools::et_eta_cut( const vector<MyParticle*> p_vector,
                     const double et_min,  const double et_max,
                     const double eta_min, const double eta_max) {

    std::vector<MyParticle*> tmp_v;
    tmp_v.clear();

    for (unsigned int i=0; i<p_vector.size(); i++){

        double et  = (p_vector[i])->tlv().Et();
        double eta = TMath::Abs((p_vector[i])->tlv().Eta());

        if(et_min>-1){
            if(et < et_min) continue; 
        }
        if(et_max>-1){
            if(et > et_max) continue; 
        }
        if(eta_min>-1){
            if(eta < eta_min) continue; 
        }
        if(eta_max>-1){
            if(eta > eta_max) continue; 
        }

        tmp_v.push_back(p_vector[i]);

    }

    return tmp_v;
}

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

//::::::::::::::::::::::::::
//:: METHOD append_vector ::
//::::::::::::::::::::::::::

std::vector<MyParticle*> 
MyTools::append_vector(std::vector<MyParticle*> first_vector, 
                       std::vector<MyParticle*> second_vector) {

    std::vector<MyParticle*> added_vector(0);
    
    for (unsigned int j=0;j<first_vector.size();j++) {
        added_vector.push_back(first_vector[j]);
    }
    
    for (unsigned int j=0;j<second_vector.size();j++) {
        added_vector.push_back(second_vector[j]);
    }
    
    return added_vector;
}

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

  //:::::::::::::::::::::::::::::::
  //:: METHOD remove_own_overlap ::
  //:::::::::::::::::::::::::::::::

  /*vector<MyParticle*> MyTools::remove_own_overlap( std::vector<MyParticle*> particle_vector, double radius ){
    double deltaR;
    double ptmax=0;
    int i_best=-1;
    vector<MyParticle*> new_vector;
    vector<MyParticle*> tmp_vector;
    new_vector.clear();
    tmp_vector.clear();
    cout <<"vector size  "<<particle_vector.size()<<endl;
    if (particle_vector.size()==1) return particle_vector;
    else {
    for (unsigned int i=1; i<particle_vector.size();i++) {
    for (unsigned int j=0; j<i;j++) {
    cout <<"i="<<i<<",  j="<<j<<endl;
    deltaR = particle_vector[i]->tlv().DeltaR(particle_vector[j]->tlv());
    if (radius>deltaR) {
    cout <<"same particles  ";
    if (particle_vector[i]->tlv().Pt()>particle_vector[j]->tlv().Pt()) {tmp_vector.push_back(particle_vector[i]);cout << i<<endl;}
    else if (particle_vector[i]->tlv().Pt()<particle_vector[j]->tlv().Pt()) {tmp_vector.push_back(particle_vector[j]);cout << j<<endl;}
    else cerr <<"WARNING: Container has particles with exactly same Pt (?)... "<<endl;
    }
    else {
    cout <<"not same particles "<<i<<"  "<<j<<endl;
    new_vector.push_back(particle_vector[i]);
    }
    }
    }
    }
    for (unsigned int i=0; i<tmp_vector.size();i++) {
    if (particle_vector[i]->tlv().Pt()>ptmax) {
    i_best=i;
    }
    }
    if (i_best!=-1) new_vector.push_back(tmp_vector[i_best]);
    return new_vector;
    }
  */


vector<MyParticle*>
MyTools::delete_doubles( vector<MyParticle*> p_vector, double delta_r){
    vector<MyParticle*> tmp_v;
    tmp_v = p_vector;
    //cout<<"delta_r: "<<delta_r<<endl;
    double dpt;
    double dr;
    dpt = 0;
    int n_doubles;
    bool double_exists;
    double_exists = false;
    n_doubles = 0;
    TLorentzVector tlv1;
    TLorentzVector tlv2;
    // loop starting from the last entry in the vector
    for (int j=tmp_v.size()-1; j>=0; j--){
        double_exists = false;
        tlv1 = (tmp_v[j])->tlv();
        // loop starting from the first entry up to the one
        // from the outer loop
        for (int k=0; k<j; k++){
            if(j==k)cerr<<"j==k in MyTools::delete_doubles"<<endl;
            tlv2 = (tmp_v[k])->tlv();
            dr = tlv1.DeltaR(tlv2);
            dpt = tlv1.Pt() - tlv2.Pt();
            if(dr<delta_r && dpt<0){
                // delete this particle
                n_doubles++;
                double_exists = true;
                //cout<<"double exists, dR: "<<dr<<", dPt: "<<dpt<<endl;
                //cout<<"pt1(to delete): "<<tlv1.Pt()<<", pt2(survives): "<<tlv2.Pt()<<endl;
            }
        }
        
        if(double_exists){
            //cout<<"removing particle at position: "<<j<<endl;
            tmp_v.erase(tmp_v.begin()+j);
        }
    }
    //cout<<"n_doubles: "<<n_doubles<<endl;
    return tmp_v; 
   
}
//*****************************************************************************

//:::::::::::::::::::::::::::::::
//:: METHOD tauRec_ETcorr      ::
//:::::::::::::::::::::::::::::::

double MyTools::tauRec_ETcorr(double et_in, double eta_in, int ntrk){
  
  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] = {15., 35., 150., 0.};
  //  if (m_gev==1) = {15000, 35000, 150000, 0.};
  double m_etapoints[m_netabins+1] = {0.25, 1.0, 2.0, 0.0};
  
  double et = et_in;
  double eta = fabs( eta_in );
  double corr = 1, fudge = 1.011;
  
  unsigned int lowpt_idx = 0;
  unsigned int loweta_idx = 0;
  
  double lowpt_frac = 0;
  double loweta_frac = 0;
  
  while ( lowpt_idx < m_nptbins-1 && et > m_ptpoints[lowpt_idx+1] ) lowpt_idx++;
  
  lowpt_frac = ( m_ptpoints[lowpt_idx+1] - et ) / ( m_ptpoints[lowpt_idx+1] -
                                                    m_ptpoints[lowpt_idx] );   
  // will be >1 if et 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;
//       std::cerr <<  et << std::endl;
//       std::cerr <<  m_nptbins     << std::endl;
//       std::cerr <<  m_ptpoints[0] << std::endl;
//       std::cerr <<  m_ptpoints[1] << std::endl;
//       std::cerr <<  m_ptpoints[2] << std::endl;
//       std::cerr <<  lowpt_idx     << std::endl;
//       std::cerr << m_ptpoints[lowpt_idx+1] << 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 (eta>2){
//       std::cerr << "---TEST---" << std::endl;
//       std::cerr <<  eta << std::endl;
//       std::cerr <<  m_netabins     << std::endl;
//       std::cerr <<  m_etapoints[0] << std::endl;
//       std::cerr <<  m_etapoints[1] << std::endl;
//       std::cerr <<  m_etapoints[2] << std::endl;
//       std::cerr <<  loweta_idx     << std::endl;
//       std::cerr << m_etapoints[loweta_idx+1] << std::endl;
//       std::cerr << m_etapoints[loweta_idx] << std::endl;
//       std::cerr <<  loweta_frac     << std::endl;
//       std::cerr << "--------" << std::endl;
//   }
  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;
//       std::cerr <<  eta << std::endl;
//       std::cerr <<  m_netabins     << std::endl;
//       std::cerr <<  m_etapoints[0] << std::endl;
//       std::cerr <<  m_etapoints[1] << std::endl;
//       std::cerr <<  m_etapoints[2] << std::endl;
//       std::cerr <<  loweta_idx     << std::endl;
//       std::cerr << m_etapoints[loweta_idx+1] << std::endl;
  }
  
  double coeff_matrix[2][2] = { {0, 0}, {0, 0} };
  
  if ( ntrk <= 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*fudge;
  
  return et;

}

//:::::::::::::::::::::::::::::::
//:: METHOD  correct_tau_et    ::
//:::::::::::::::::::::::::::::::

vector<MyParticle*> MyTools::correct_tau_et( std::vector<MyParticle*> tau_vector){

        vector<MyParticle*> new_vector;
        new_vector.clear();
        MyTauJet *tmp_tau = new MyTauJet;
        for (unsigned int k=0; k<tau_vector.size();k++) {

          double tau_et=((MyTauJet*)(tau_vector[k]))->ethadcalib()
            +((MyTauJet*)(tau_vector[k]))->etemcalib();
          double tau_eta=(tau_vector[k])->tlv().Eta();
          double tau_phi=(tau_vector[k])->tlv().Phi();
          int tau_ntrk=((MyTauJet*)(tau_vector[k]))->nb_tracks();
          
          
          double tau_et_recalib=tauRec_ETcorr(tau_et, tau_eta, tau_ntrk);

          double tau_e_recalib=tau_et_recalib*cosh(tau_eta);
          double tau_p2_recalib=tau_e_recalib*tau_e_recalib-1.77699*1.77699;
          double tau_p_recalib=sqrt(tau_p2_recalib);
          double tau_pt_recalib=tau_p_recalib/cosh(tau_eta);;

          TLorentzVector new_tlv;
          new_tlv.SetPtEtaPhiE(tau_pt_recalib,tau_eta,tau_phi,tau_e_recalib);

          tmp_tau=((MyTauJet*)tau_vector[k])->Clone();
          tmp_tau->set_tlv(new_tlv);
          new_vector.push_back(tmp_tau);

        }
        return new_vector;

}

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

//::::::::::::::::::::::::
//:: METHOD PrintVector ::
//::::::::::::::::::::::::

void
MyTools::PrintVector( vector<MyParticle*> p_vector, std::string option){
    
    for (unsigned int k=0; k<p_vector.size(); k++) {
        (p_vector[k])->Print(option);
    } 

}



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

//::::::::::::::::::::::::::::::::::
//:: METHOD particles_from_vertex ::
//::::::::::::::::::::::::::::::::::

vector<MyParticle*> 
MyTools::particles_from_vertex(const vector<MyParticle*> & p_vector, const int vertex) {
    
    vector<MyParticle*> particles;

    for (unsigned int k=0; k<p_vector.size(); k++) {
        
        if(p_vector[k]->vertexIndex() == vertex){
            particles.push_back(p_vector[k]);
        }
    }
    
    return particles;
    
}

//:::::::::::::::::::::::::::::::
//:: METHOD change_btag_weight ::
//:::::::::::::::::::::::::::::::

double MyTools::change_btag_weight(double current_btag_weight, double delta_weight){

  double new_weight;
  new_weight = current_btag_weight + delta_weight;

  return new_weight;

}

//:::::::::::::::::::::::::
//:: METHOD scale_b_tags ::
//:::::::::::::::::::::::::

void MyTools::scale_b_tags(std::vector<MyParticle*> * jets, int flavour, double delta_bweight){

  for (unsigned int il=0;il<jets->size();il++){

    double bweight=((MyParticleJet*)((*jets)[il]))->b_tag_weight();
    int pdgId=((MyParticleJet*)((*jets)[il]))->truth_flavour();
    double new_bweight=bweight;

    if (pdgId==flavour){
      new_bweight=change_btag_weight(bweight,delta_bweight);
    }

    ((MyParticleJet*)((*jets)[il]))->set_b_tag_weight(new_bweight);

//     cout<<"flavour="<<pdgId<<" old="<<bweight<<" delta="<<delta_bweight<<" new="<<new_bweight
//      <<" new_pack="<<((MyParticleJet*)((*jets)[il]))->b_tag_weight()<<"\n";

  }

}

---