TMBJet.cpp

Go to the documentation of this file.

#include "tmb_tree/TMBJet.hpp"
#include <iostream>
#include "TRandom.h" //for btag check
#include "TTree.h"
#include "TBranch.h"

#include <cassert>
#include <stdexcept>

ClassImp(TMBJet);

TMBJet::TMBJet()
{
}

TMBJet::TMBJet      (Float_t E, Float_t pT,   Float_t phi,  Float_t eta)
{
    Set0   (E, pT, phi, eta);
}

TMBJet::~TMBJet()
{
}

void TMBJet::Set0    (Float_t E, Float_t pT,   Float_t phi,  Float_t eta)
{

    SetPtEtaPhiE(pT, eta, phi, E);

    //  alwyas initialize corrections to 0 
    _jes_C=0.0;
    _jesMU_C=0.0;
    _jes_dC_stat=0.0;
    _jesMU_dC_stat=0.0;
    _jes_dC_stat_up=0.0;
    _jesMU_dC_stat_up=0.0;
    _jes_dC_stat_down=0.0;
    _jesMU_dC_stat_down=0.0;
    _jes_dC_sys=0.0;
    _jesMU_dC_sys=0.0;
    _jes_dC_sys_up=0.0;
    _jesMU_dC_sys_up=0.0;
    _jes_dC_sys_down=0.0;
    _jesMU_dC_sys_down=0.0;
    _jes_metC=0.0;
    _jesMU_metC=0.0;
    _jes_metdC_stat=0.0;
    _jesMU_metdC_stat=0.0;
    _jes_metdC_sys=0.0;
    _jesMU_metdC_sys=0.0;
    _smear_coeff = 1.0;
    _smear_coeffMU = 1.0;
    _isGood = 1;
    _isL1Conf = 0;
    _isBad = 0;
    _isNoise = 0;
    _isEM = 0;
    _hasMU = 0;
    _isSmeared = 0;
    _isSmearedMU = 0;
    _actas = 0;

    ClearAllBTags();
}

void TMBJet::ActAsUnCorrected() const
{
  if (_actas==kUnCorrected) return;
  const_cast<TMBJet&>(*this).SetPtEtaPhiE(_uncorrected.Pt(), _uncorrected.Eta(), _uncorrected.Phi(), _uncorrected.E());
  _actas = kUnCorrected;
  return;
}

void TMBJet::ActAsJESCorrected() const
{
  if (_actas==kJESCorrected) return;
  const_cast<TMBJet&>(*this).SetPEtaPhiE(_uncorrected.P()*_jes_C, _uncorrected.Eta(), _uncorrected.Phi(), _uncorrected.E()*_jes_C);
  _actas = kJESCorrected;
  return;
}

void TMBJet::ActAsJESCorrectedShiftedPlus() const
{
  if (_actas==kJESCorrectedShiftedPlus) return;
  const_cast<TMBJet&>(*this).SetPEtaPhiE(_uncorrected.P()*(_jes_C+sqrt(_jes_dC_sys_up*_jes_dC_sys_up+_jes_dC_stat_up*_jes_dC_stat_up)),
                 _uncorrected.Eta(), _uncorrected.Phi(), 
                 _uncorrected.E()*(_jes_C+sqrt(_jes_dC_sys_up*_jes_dC_sys_up+_jes_dC_stat_up*_jes_dC_stat_up)));
  _actas = kJESCorrectedShiftedPlus;
  return;
}

void TMBJet::ActAsJESCorrectedShiftedMinus() const
{
  if (_actas==kJESCorrectedShiftedMinus) return;
  const_cast<TMBJet&>(*this).SetPEtaPhiE(_uncorrected.P()*(_jes_C-sqrt(_jes_dC_sys_down*_jes_dC_sys_down+_jes_dC_stat_down*_jes_dC_stat_down)), 
                _uncorrected.Eta(), _uncorrected.Phi(), 
                _uncorrected.E()*(_jes_C-sqrt(_jes_dC_sys_down*_jes_dC_sys_down+_jes_dC_stat_down*_jes_dC_stat_down)));
  _actas = kJESCorrectedShiftedMinus;
  return;
}

void TMBJet::ActAsJESMUCorrected() const
{
  if (_actas==kJESMUCorrected) return;
  const_cast<TMBJet&>(*this).SetPEtaPhiE(_uncorrected.P()*_jesMU_C, _uncorrected.Eta(), _uncorrected.Phi(), _uncorrected.E()*_jesMU_C);
  _actas = kJESMUCorrected;
  return;
}

void TMBJet::ActAsJESMUCorrectedShiftedPlus() const
{
  if (_actas==kJESMUCorrectedShiftedPlus) return;
  const_cast<TMBJet&>(*this).SetPEtaPhiE(_uncorrected.P()*(_jesMU_C+sqrt(_jesMU_dC_sys_up*_jesMU_dC_sys_up+_jesMU_dC_stat_up*_jesMU_dC_stat_up)),
                                         _uncorrected.Eta(), _uncorrected.Phi(), 
                                         _uncorrected.E()*(_jesMU_C+sqrt(_jesMU_dC_sys_up*_jesMU_dC_sys_up+_jesMU_dC_stat_up*_jesMU_dC_stat_up)));
  _actas = kJESMUCorrectedShiftedPlus;
  return;
}

void TMBJet::ActAsJESMUCorrectedShiftedMinus() const
{
  if (_actas==kJESMUCorrectedShiftedMinus) return;
  const_cast<TMBJet&>(*this).SetPEtaPhiE(_uncorrected.P()*(_jesMU_C-sqrt(_jesMU_dC_sys_down*_jesMU_dC_sys_down+_jesMU_dC_stat_down*_jesMU_dC_stat_down)), 
                                         _uncorrected.Eta(), _uncorrected.Phi(), 
                                         _uncorrected.E()*(_jesMU_C-sqrt(_jesMU_dC_sys_down*_jesMU_dC_sys_down+_jesMU_dC_stat_down*_jesMU_dC_stat_down)));
  _actas = kJESMUCorrectedShiftedMinus;
  return;
}

void TMBJet::ActAsSmeared() const
{
  if (_actas==kSmeared) return;
  const_cast<TMBJet&>(*this).SetPEtaPhiE(_uncorrected.P()*_jes_C*_smear_coeff, _uncorrected.Eta(), _uncorrected.Phi(), _uncorrected.E()*_jes_C*_smear_coeff);
  _actas = kSmeared;
  return;
}

void TMBJet::ActAsSmearedMU() const
{
  if (_actas==kSmearedMU) return;
  const_cast<TMBJet&>(*this).SetPEtaPhiE(_uncorrected.P()*_jesMU_C*_smear_coeffMU, _uncorrected.Eta(), _uncorrected.Phi(), _uncorrected.E()*_jesMU_C*_smear_coeffMU);
  _actas = kSmearedMU;
  return;
}

void TMBJet::SetActAs(const Int_t ActAs) const
{
  if (ActAs == TMBJet::kUnCorrected) {
    ActAsUnCorrected();
  } else if (ActAs == TMBJet::kJESCorrected) {
    ActAsJESCorrected();
  } else if (ActAs == TMBJet::kJESCorrectedShiftedPlus) {
    ActAsJESCorrectedShiftedPlus();
  } else if (ActAs == TMBJet::kJESCorrectedShiftedMinus) {
    ActAsJESCorrectedShiftedMinus();
  } else if (ActAs == TMBJet::kJESMUCorrected) {
    ActAsJESMUCorrected();
  } else if (ActAs == TMBJet::kJESMUCorrectedShiftedPlus) {
    ActAsJESMUCorrectedShiftedPlus();
  } else if (ActAs == TMBJet::kJESMUCorrectedShiftedMinus) {
    ActAsJESMUCorrectedShiftedMinus();
  } else if (ActAs == TMBJet::kSmeared) {
    ActAsSmeared();
  } else if (ActAs == TMBJet::kSmearedMU) {
    ActAsSmearedMU();
  } else {
    throw std::runtime_error( Form ("TMBJet: Can't reset jet's ActAs to unknown JES type (%d)", ActAs));
  }
  return;
}

void TMBJet::Set1    (Float_t q, Float_t dphi, Float_t deta, 
                      Float_t emf,  Float_t em1f, Float_t em2f, Float_t em3f,
                      Float_t ccmg, Float_t icdf, Float_t ecmg, Float_t icrf,
                      Float_t fh1f, Float_t fh2f, Float_t fh3f, Float_t chf,
                      Float_t emcc, Float_t hadcc,Float_t emec, Float_t hadec,
                      Float_t hot,  Float_t mxET,
                      Float_t cpsE, Float_t etaW, Float_t phiW, 
                      Float_t sET,  Float_t vPT,  Float_t iPT,  Float_t seedET,
                      Float_t pxCH, Float_t pyCH, Float_t pzCH)
{
    _q     = q;
    _dphi  = dphi; 
    _deta  = deta; 
    _emf   = emf; 
    _em1f  = em1f; 
    _em2f  = em2f; 
    _em3f  = em3f;
    _ccmg  = ccmg;
    _icdf  = icdf; 
    _ecmg  = ecmg;
    _icrf  = icrf; 
    _fh1f  = fh1f; 
    _fh2f  = fh2f; 
    _fh3f  = fh3f; 
    _emcc  = emcc;
    _hadcc = hadcc;
    _emec  = emec;
    _hadec = hadec;
    _chf   = chf; 
    _hot   = hot; 
    _mxET  = mxET; 
    _cpsE  = cpsE; 
    _etaW  = etaW; 
    _phiW  = phiW; 
    _sET   = sET; 
    _vPT   = vPT; 
    _iPT   = iPT; 
    _seedET= seedET; 
    _pxCH = pxCH;
    _pyCH = pyCH;
    _pzCH = pzCH;
}


void TMBJet::Set2    (Int_t ntrk, Int_t nps,  Int_t Nitems, 
                      Int_t n90,  Int_t split_merge_word, 
                      const char *algoname, const TRefArray *tracks, TRef vtxref)
{
    _ntrk  = ntrk; 
    _nps   = nps; 
    _Nitems= Nitems; 
    _n90   = n90;
    _split_merge_word =  split_merge_word;
    _algoname = algoname;              
    _vtxref   = vtxref;

    _tracks.Clear();
    if (tracks) {
        for (Int_t i=0; i<tracks->GetLast()+1; i++) {
            _tracks.Add( tracks->At(i) );
        }
    }
}

void TMBJet::SetJetShapes(const std::vector<std::pair<Char_t, UChar_t> >& jetTowers, 
                          const std::vector<Float_t>& jetShapes)
{
  _jetTowers=jetTowers;
  _jetShapes=jetShapes;
}

void TMBJet::SetL1(Float_t l1set, Float_t l1pt, Float_t l1em)
{
    _l1set = l1set;
    _l1pt  = l1pt;
    _l1em  = l1em;
}

void TMBJet::SetCorr (Float_t jes_C, Float_t jes_dC_stat, Float_t jes_dC_sys,
                      Float_t jes_metC, Float_t jes_metdC_stat, Float_t jes_metdC_sys)
{
    _jes_C=jes_C;
    _jes_dC_stat=jes_dC_stat;
    _jes_dC_sys=jes_dC_sys;
    _jes_metC=jes_metC;
    _jes_metdC_stat=jes_metdC_stat;
    _jes_metdC_sys=jes_metdC_sys;

 
    switch(_actas) {
    case kJESCorrected:
        _actas = kUnCorrected;
        ActAsJESCorrected();
        break;
    case kJESCorrectedShiftedPlus:
        _actas = kUnCorrected;
        ActAsJESCorrectedShiftedPlus();
        break;
    case kJESCorrectedShiftedMinus:
        _actas = kUnCorrected;
        ActAsJESCorrectedShiftedMinus();
        break;
    default:
        break;
    }    
}

void TMBJet::SetCorr (Float_t jes_C, Float_t jes_dC_stat_up, Float_t jes_dC_sys_up, Float_t jes_dC_stat_down, Float_t jes_dC_sys_down,
                      Float_t jes_metC, Float_t jes_metdC_stat, Float_t jes_metdC_sys)
{
    _jes_C=jes_C;
    _jes_dC_stat_down=jes_dC_stat_down;
    _jes_dC_sys_down=jes_dC_sys_down;
    _jes_dC_stat_up=jes_dC_stat_up;
    _jes_dC_sys_up=jes_dC_sys_up;
    _jes_dC_stat = _jes_dC_stat_up >= _jes_dC_stat_down ? _jes_dC_stat_up : _jes_dC_stat_down ;
    _jes_dC_sys  = _jes_dC_sys_up >= _jes_dC_sys_down ? _jes_dC_sys_up : _jes_dC_sys_down ;
    _jes_metC=jes_metC;
    _jes_metdC_stat=jes_metdC_stat;
    _jes_metdC_sys=jes_metdC_sys;
 
    switch(_actas) {
    case kJESCorrected:
        _actas = kUnCorrected;
        ActAsJESCorrected();
        break;
    case kJESCorrectedShiftedPlus:
        _actas = kUnCorrected;
        ActAsJESCorrectedShiftedPlus();
        break;
    case kJESCorrectedShiftedMinus:
        _actas = kUnCorrected;
        ActAsJESCorrectedShiftedMinus();
        break;
    default:
        break;
    }

}

void TMBJet::SetCorrMU(Float_t jes_C, Float_t jes_dC_stat, Float_t jes_dC_sys,
                      Float_t jes_metC, Float_t jes_metdC_stat, Float_t jes_metdC_sys)
{
    _jesMU_C=jes_C;
    _jesMU_dC_stat=jes_dC_stat;
    _jesMU_dC_sys=jes_dC_sys;
    _jesMU_metC=jes_metC;
    _jesMU_metdC_stat=jes_metdC_stat;
    _jesMU_metdC_sys=jes_metdC_sys;


    switch(_actas) {
    case kJESMUCorrected:
        _actas = kUnCorrected;
        ActAsJESMUCorrected();
        break;
    case kJESMUCorrectedShiftedPlus:
        _actas = kUnCorrected;
        ActAsJESMUCorrectedShiftedPlus();
        break;
    case kJESMUCorrectedShiftedMinus:
        _actas = kUnCorrected;
        ActAsJESMUCorrectedShiftedMinus();
        break;
    default:
        break;
    }   
}

void TMBJet::SetCorrMU(Float_t jes_C, Float_t jes_dC_stat_up, Float_t jes_dC_sys_up, Float_t jes_dC_stat_down, Float_t jes_dC_sys_down,
                      Float_t jes_metC, Float_t jes_metdC_stat, Float_t jes_metdC_sys)
{
    _jesMU_C=jes_C;
    _jesMU_dC_stat_up=jes_dC_stat_up;
    _jesMU_dC_sys_up=jes_dC_sys_up;
    _jesMU_dC_stat_down=jes_dC_stat_down;
    _jesMU_dC_sys_down=jes_dC_sys_down;
    _jesMU_dC_stat = _jesMU_dC_stat_up >= _jesMU_dC_stat_down ? _jesMU_dC_stat_up : _jesMU_dC_stat_down ;
    _jesMU_dC_sys  = _jesMU_dC_sys_up >= _jesMU_dC_sys_down ? _jesMU_dC_sys_up : _jesMU_dC_sys_down ;
    _jesMU_metC=jes_metC;
    _jesMU_metdC_stat=jes_metdC_stat;
    _jesMU_metdC_sys=jes_metdC_sys;


    switch(_actas) {
    case kJESMUCorrected:
        _actas = kUnCorrected;
        ActAsJESMUCorrected();
        break;
    case kJESMUCorrectedShiftedPlus:
        _actas = kUnCorrected;
        ActAsJESMUCorrectedShiftedPlus();
        break;
    case kJESMUCorrectedShiftedMinus:
        _actas = kUnCorrected;
        ActAsJESMUCorrectedShiftedMinus();
        break;
    default:
        break;
    }
}

void TMBJet::SetFlag(Int_t isGood, Int_t isL1Conf, Int_t isBad, Int_t isNoise, Int_t isEM, Int_t hasMU, Int_t isSmeared, Int_t isSmearedMU)
{
  _isGood      = isGood;
  _isL1Conf    = isL1Conf;
  _isBad       = isBad;
  _isNoise     = isNoise;
  _isEM        = isEM;
  _hasMU       = hasMU;
  _isSmeared   = isSmeared;
  _isSmearedMU = isSmearedMU;

}
                                                                                                                           
void TMBJet::SetUnCorr(Float_t p, Float_t eta, Float_t phi, Float_t E)
{
  TMBLorentzVector temp(p,eta,phi,E,kPEtaPhiE);
  _uncorrected = temp;
}

void TMBJet::SetSmearCoeff(Float_t coeff, Float_t coeffMU)
{
    _smear_coeff = coeff;
    _smear_coeffMU = coeffMU;

    switch(_actas) {
    case kSmeared:
        _actas = kUnCorrected;
        ActAsSmeared();
        break;
    case kSmearedMU:
        _actas = kUnCorrected;
        ActAsSmearedMU();
        break;
    default:
        break;
    }
}

void TMBJet::SetJetID(Float_t cpf0, Float_t cpf, Int_t ntrkMultiplicity0,
                      Int_t ntrkMultiplicity, TRef cpfvtxref)
{
    _cpf0 = cpf0;
    _cpf  = cpf;
    _ntrkMultiplicity0 = ntrkMultiplicity0;
    _ntrkMultiplicity  = ntrkMultiplicity;
    _cpfvtxref = cpfvtxref;
}

void TMBJet::SetCellList(TMBCellContainer *cont, 
                         const std::vector<TMBCellContainer::Index_t>& cells)
{
    _container = cont;
    _cells = cells;
}


UInt_t TMBJet::ncells() const
{
    return _cells.size();
}

const TMBCaloCell *TMBJet::GetCaloCell(UInt_t icell) const
{

    assert(icell < (UInt_t)_cells.size());
    if(TMBCellContainer *c = _container.IsValid() ? (TMBCellContainer *)_container.GetObject() : 0) {
        assert(_cells[icell] < c->NumCells());
        return c->GetCell(_cells[icell]);
    } else {
        return 0;
    }    
}

/*
namespace {

    //Read the BTag branch being refered to
    bool ReadBTagBranch(const TString& algo, const TString& bname){
        TTree *tree = (TTree*)gROOT->FindObject("TMBTree");
        if (tree){
            Long64_t current = tree->GetReadEntry();
            //std::cout<<"entry="<<current<<std::endl;
            TBranch *b =tree->GetBranch("BTag_corr"+algo+"_"+bname);
            if (current>=0 && b) {
                Long64_t other_current = b->GetReadEntry();
                //std::cout<<"other entry="<<other_current<<std::endl;
                if (other_current!=current) {
                    //std::cout<<"reading entry "<<current<<" from branch "<<"BTag_corrJCCB_"+i->first<<std::endl;
                    b->GetEntry(current);
                }
            }
            return true;
        }
        return false;
    }
}
*/

// Return a pointer the TMBBTag object associated with this jet 
// for a particular tag type ("SVT", etc.) and cut ("LOOSE", etc.). 
// Note, you may have to do a dynamic_cast to access the subclass info, as for TMBBTagSVT objects, 
// for instance:
//    TMBBTagSVT* tag = dynamic_cast<TMBBTagSVT*> GetBTag("NN","MEDIUM");
// Also, you should check that the pointer returned is not NULL! 
// It will be if no tag info for that type of tag is associated with the jet.
const TMBBTag* TMBJet::GetBTag(const TString &type, const TString &cut) const {
    TString typecut(type+"_"+cut);
    return GetBTag(typecut);
}

const TMBBTag* TMBJet::GetBTag(const TString &typecut) const {

    //first get the btag check integer
    int check = -1;
    std::map< TString, int >::const_iterator c = _btagdictcheck.find(typecut);
    if (_btagdictcheck.end()!=c) {
        check = c->second;
    }

    std::map< TString, TRef >::const_iterator i = _btagdict.find(typecut);
    if (_btagdict.end()!=i && i->second.IsValid()) {
    
        //First make sure we read the BTag branch being referenced
        //ReadBTagBranch(_algoname, i->first);
    
        const TMBBTag* tag = (const TMBBTag*)(i->second.GetObject());
        if (tag) {
          if (tag->GetId()!=check){
            std::cout<<"WARNING! ID check failed in GetBTag!"<<std::endl;
          }
        }
        return tag;
    }
    else {
        return 0;
    }
}

const TMBBTag* TMBJet::GetBTag(const char *type, const char *cut) const
{
    return GetBTag (TString (type), TString (cut));
}

//Get the first BTag object available 
//(useful when you don't care about the algorithm itself, like for taggability or mc_flavor...
const TMBBTag* TMBJet::GetFirstBTag() const 
{  
    std::map< TString, TRef >::const_iterator i = _btagdict.begin();
    while (i != _btagdict.end()) {
      const TMBBTag *result = GetBTag(i->first);
      if (result != 0) {
        return result;
      }
      i++;
    }
    return 0;
}

//Print the list of associated, valid BTags in the map
bool TMBJet::btag_print() const 
{
    std::cout<<"btags for this jet:";
    bool ret=true;
    std::map< TString, TRef >::const_iterator i = _btagdict.begin();
    while (_btagdict.end()!=i) {

      //First make sure we read the BTag branch being referenced
      //ReadBTagBranch(_algoname, i->first);
      
      std::cout<<" "<<i->first;
      if (!i->second.IsValid()) {std::cout<<"-invalid! "; ret=false;}
      
      const TMBBTag* tag = (const TMBBTag*)(i->second.GetObject());
      if (!tag) {std::cout<<" <-tag==0!"; ret=false;}
      else {
        
        //check btag int
        std::cout<<" "<<tag->GetId()<<"="<<btagdictcheck(i->first);
        if (tag->GetId()!=btagdictcheck(i->first)) {std::cout<<" <-mismatch!"; ret=false;}      
        
        //check jet ref
        const TMBJet* jet = tag->GetJet();
        if (!jet) std::cout<<" bad jet ref";
        else {
          
          //do the Pt's match?
          std::cout<<" "<<jet->Pt()<<"="<<this->Pt();
          if (fabs(jet->Pt()-this->Pt())>0.000001) {std::cout<<" <-mismatch!"; ret=false;}
          
          //check btag int
          std::cout<<" "<<tag->GetId()<<"="<<jet->btagdictcheck(i->first);
          if (tag->GetId()!=jet->btagdictcheck(i->first)) {std::cout<<" <-mismatch!"; ret=false;}
          
        }
        
      }
      ++i;
    }
    std::cout<<std::endl;
    return ret;
}

void TMBJet::AddBTag(const TString &type, const TString &cut, TMBBTag* btag)
{
    TString typecut(type+"_"+cut);

    //This is a check.... generate a random integer, and store it in both the btag and the map of the jet.
    //Later you can check that they still match!
    int r = gRandom->Integer(1000000000);
    btag->SetId(r);
    _btagdictcheck[typecut] = r;

    //This is the real map entry
    _btagdict[typecut] = btag;
}

void TMBJet::ClearAllBTags()
{
    _btagdict.clear();
    _btagdictcheck.clear();
}

void TMBJet::ClearBTag (const TString &type, const TString &cut)
{
    std::map<TString, TRef>::iterator itr = _btagdict.find(type+"_"+cut);
    if (itr != _btagdict.end()) {
        _btagdict.erase(itr);
        //std::cout<<"Erasing tag for "<<type<<" _ "<<cut<<" - size = "<<_btagdict.size()<<std::endl;
    }
    std::map<TString, int>::iterator itr2 = _btagdictcheck.find(type+"_"+cut);
    if (itr2 != _btagdictcheck.end()) {
        _btagdictcheck.erase(itr2);
        //std::cout<<"Erasing check tag for "<<type<<" _ "<<cut<<" - size = "<<_btagdictcheck.size()<<std::endl;
    }
}

int TMBJet::btagdictcheck(const char* t) const 
{
    if (_btagdictcheck.size()==0) return -2;
    TString s(t);

    std::map< TString, int >::const_iterator i = _btagdictcheck.find(s);
    if (i!=_btagdictcheck.end()) return i->second;
  
    //now look for a close match... something ending in the TString key, like the algo name
    for (i=_btagdictcheck.begin(); i!=_btagdictcheck.end(); ++i){
        if (s.EndsWith(i->first)) return i->second;
    }

    // If there is only one associated btag object, return that.
    if (_btagdictcheck.size() == 1)
      return _btagdictcheck.begin()->second;
  
    //ok, nothing's close
    return -1;
}

int TMBJet::taggable() const 
{ 
    const TMBBTag* svt = GetFirstBTag();
    if (svt) return svt->is_taggable();
    else return -1;
}

int TMBJet::mc_flavor() const 
{ 
    const TMBBTag* svt = GetFirstBTag();
    if (svt) return svt->mc_flavor();
    else return -1;
}

---