MultivariateAnalysis/src/LJetsTopoVars.cc

#include "LJMet/MultivariateAnalysis/interface/LJetsTopoVars.h"
#include "LJMet/MultivariateAnalysis/interface/TopTopologicalVariables.h"
#include "LJMet/MultivariateAnalysis/interface/AnglesUtil.h"

#include "TMatrixDSymEigen.h"

#include <stdexcept>

#include "DataFormats/PatCandidates/interface/Jet.h"
#include "DataFormats/PatCandidates/interface/Electron.h"
#include "DataFormats/PatCandidates/interface/Muon.h"
#include "DataFormats/PatCandidates/interface/MET.h"

using namespace std;
//using namespace cafe;
using namespace top_cafe;


namespace
{
   bool moreThan(const pat::Jet& a, const pat::Jet& b)
   {
      return a.pt() > b.pt();    
   }
}


int LJetsTopoVars::setEvent(const edm::Event& event, double min_dr_jet_lepton)
{
  using namespace edm;

  int removed_jets = 0;

  Handle< vector< pat::Jet > > jets;
  Handle< vector< pat::MET > > met;
  Handle< vector< pat::Electron > > electrons;
  Handle< vector< pat::Muon > > muons;

  event . getByLabel( m_jetBranch, jets );
  event . getByLabel( m_metBranch, met );
  if (m_isMuon) event . getByLabel( m_leptonBranch, muons );
  else event . getByLabel( m_leptonBranch, electrons );

  //cout << endl << "=====> LJetsTopoVars::setEvent(): number of jets = " << jets -> size() << endl;
  
  m_jets.clear();
  
  eigenval.ResizeTo(3);
  eigenval.Zero();
  

  if (met->size()>0){
    m_met = TMBLorentzVector(met->begin()->pt(),met->begin()->eta(),met->begin()->phi(),met->begin()->energy(),TMBLorentzVector::kPtEtaPhiE);
  }
  else{
    cout << "LJetsTopoVars::setEvent(): no MET in this event!" << endl;
  }
   
  if(m_isMuon){
    if (muons->size()>0){
      m_lepton = TMBLorentzVector(muons->begin()->pt(),muons->begin()->eta(),muons->begin()->phi(),muons->begin()->energy(),TMBLorentzVector::kPtEtaPhiE);
      //_muon = &(event.getCollection<TMBMuon>(m_leptonBranch.c_str())[0]);
    }
  }
  else{
    if (electrons->size()>0){
      m_lepton = TMBLorentzVector(electrons->begin()->pt(),electrons->begin()->eta(),electrons->begin()->phi(),electrons->begin()->energy(),TMBLorentzVector::kPtEtaPhiE);
      //_electron = &(event.getCollection<TMBEMCluster>(m_leptonBranch.c_str())[0]);
    }
  }


  // loop over first four jets (that are not too close to the lepton!!!)
  //vector<pat::Jet>::const_iterator jet;
  //int nMax = jets->size() > 4 ? 4 : jets->size();
  //jet = jets->begin();
  //for (int i=0; i<nMax; ++i){
  for (vector<pat::Jet>::const_iterator jet=jets->begin(); (jet!=jets->end()) && (m_jets.size()!=4); jet++){
    //cout << "LJetsTopoVars::setEvent(): jet pt() = " << jet -> pt() << endl;
    TMBLorentzVector _j(jet->pt(),jet->eta(),jet->phi(),jet->energy(),TMBLorentzVector::kPtEtaPhiE);
    if (m_lepton.DeltaR(_j) > min_dr_jet_lepton){
      m_jets.push_back(_j);
      //jet++;
    }
    else{
      removed_jets++;
    }
  }

        //original method had as given:
        //        std::vector<TheJetClass*> selectedJets, 
        //    TLorentzVector selectedLepton,
        //    double nu_px (MET px OK?)
        //    double nu_py (MET py OK?)
        double nu_px = m_met.Px();
        double nu_py = m_met.Py();

        //set all OK flags to FALSE;
        _htOK = false;
        _evtTopoOK = false;
        _ktOK = false;
        _mtOK = false;

        //
        // calculate neutrino lorentz vector (from Tobi's TopSvtAnalysis)
        //
        double nu_pz = 0.;
        double nu_e  = sqrt(pow(nu_px,2)+pow(nu_py,2));

        double Mw    = 80.43;  // NGO fix this!(read from one place)
        double l_px  = m_lepton.Px();
        double l_py  = m_lepton.Py();
        double l_pz  = m_lepton.Pz();
        double l_pt  = m_lepton.Pt();
        double l_e   = m_lepton.E();
        double Mt    = sqrt(pow(l_pt+nu_e ,2)-
                       pow(l_px+nu_px,2)-
                       pow(l_py+nu_py,2));

        double A;
        if (Mt<Mw) A = pow(Mw,2)/2.;
        else       {           // assume Mt=Mw, rescale MET accordingly (NGO???)
                A = pow(Mt,2)/2.;
                double k = nu_e*l_pt - nu_px*l_px - nu_py*l_py;
                k = (k == 0. ? 0.00001 : k);
                double scf = 0.5*pow(Mw,2)/k ;
                nu_px *= scf;
                nu_py *= scf;
                nu_e = sqrt(pow(nu_px,2)+pow(nu_py,2));
        }    

        double B = nu_px*l_px + nu_py*l_py;
        double C = TMath::Max(1. + pow(nu_e,2) * (pow(l_pz,2)-pow(l_e,2)) / pow(A+B,2)  , 0.);
        C = sqrt(C);
        double S1= (-(A+B)*l_pz + (A+B)*l_e*C) / (pow(l_pz,2)-pow(l_e,2));
        double S2= (-(A+B)*l_pz - (A+B)*l_e*C) / (pow(l_pz,2)-pow(l_e,2));

        // choose solution with smallest |l_pz| a la Run I
        nu_pz = fabs (S1) < fabs (S2) ? S1 : S2 ;

        //NGO: NOTE: neutrino PX, PY are not necessarily metPX, metPY any more!!!
        _neutrino.SetPxPyPzE(nu_px,nu_py,nu_pz,nu_e);

  /*****
  *****/
        return removed_jets;
}

double LJetsTopoVars::aplanarity() const
{
   vector<TMBLorentzVector> objects(m_jets);
   objects.push_back(m_lepton);
   TopTopologicalVariables jetsPlusLepton(objects);
   return jetsPlusLepton.Aplanarity();
}

double LJetsTopoVars::centrality() const
{
   TopTopologicalVariables jets(m_jets);
   return jets.Centrality();
}

double LJetsTopoVars::sphericity() const
{
   vector<TMBLorentzVector> objects(m_jets);
   objects.push_back(m_lepton);
   TopTopologicalVariables jetsPlusLepton(objects);
   return jetsPlusLepton.Sphericity();
}

double LJetsTopoVars::ht() const
{
   TopTopologicalVariables jets(m_jets);
   return jets.Ht();
}

double LJetsTopoVars::htpluslepton() const
{
   vector<TMBLorentzVector> objects(m_jets);
   objects.push_back(m_lepton);
   TopTopologicalVariables jetsPlusLepton(objects);
   return jetsPlusLepton.Ht();
}

double LJetsTopoVars::methtpluslepton() const
{
   vector<TMBLorentzVector> objects(m_jets);
   objects.push_back(m_lepton);
   objects.push_back(m_met);
   TopTopologicalVariables metjetsPlusLepton(objects);
   return metjetsPlusLepton.Ht();
}

double LJetsTopoVars::h() const
{
   TopTopologicalVariables jets(m_jets);
   return jets.H();
}


double LJetsTopoVars::ktMinPrime() const
{
   TopTopologicalVariables jets(m_jets);
   float ktmin = jets.KtMin();
   float etw = m_met.Pt() + m_lepton.Pt();
   return ktmin/etw;
}

double LJetsTopoVars::dphiLepMet() const
{
   return kinem::delta_phi(m_met.Phi(), m_lepton.Phi());
}

double LJetsTopoVars::minDijetMass() const
{
   TopTopologicalVariables jets(m_jets);
   return jets.MinimumPairMass();
}

double LJetsTopoVars::maxJetEta() const 
{
        double jetEta = 0;
        for (unsigned int i=0; i<m_jets.size(); i++) {
          if(TMath::Abs(m_jets.at(i).Eta()) > TMath::Abs(jetEta) ) jetEta = TMath::Abs(m_jets.at(i).Eta());
        }
        return jetEta;
}


double LJetsTopoVars::Et3() const 
{
        double Et3 = 0;
        for (unsigned int i=2; i<m_jets.size(); i++) {
          Et3+=m_jets.at(i).Pt();
        }
        return Et3;
}

double LJetsTopoVars::minDijetDeltaR() const
{

  int nJet = m_jets.size();

  double dRmin = 9999.;
  double eTmin = 9999.;
  for(int i=0;i<nJet-1;i++){
    for(int j=i+1;j<nJet;j++){
      double dR = m_jets[i].DeltaR(m_jets[j]);
      if(dR<dRmin){
        dRmin = dR;
        eTmin = std::min(m_jets[i].Pt(),m_jets[j].Pt());
      }
    }
  }
  if(dRmin>100.) {dRmin=0.;}
  
  return dRmin;
}


double LJetsTopoVars::Hz() {
        vector<TMBLorentzVector> objects;
        objects.assign(m_jets.begin(), m_jets.end());
        objects.push_back(m_lepton);
        objects.push_back(_neutrino);
        double pz = 0;
        for (vector<TMBLorentzVector>::iterator obj = objects.begin(); obj!=objects.end(); ++obj) pz += abs((*obj).Pz());
        return pz;
}

double LJetsTopoVars::HT2() {
        vector<TMBLorentzVector> objects;
        objects.assign(++m_jets.begin(), m_jets.end());
        TopTopologicalVariables topo(objects);
        return topo.Ht();
}

double LJetsTopoVars::HT2prime() {
        return HT2()/Hz();
}

double LJetsTopoVars::W_MT() {
        vector<TMBLorentzVector> objects;
        //objects.push_back(_neutrino);  //_neutrino was made with W mass constraint; use MET instead
        objects.push_back(m_met);
        objects.push_back(m_lepton);
        TopTopologicalVariables topo(objects);
        return topo.TransverseMass();
}

double LJetsTopoVars::W_Pt() {
        vector<TMBLorentzVector> objects;
        //objects.push_back(_neutrino);  //_neutrino was made with W mass constraint; use MET instead
        objects.push_back(m_met);
        objects.push_back(m_lepton);
        TopTopologicalVariables topo(objects);
        return topo.Pt();
}

double LJetsTopoVars::W_M() {
        vector<TMBLorentzVector> objects;
        objects.push_back(_neutrino);  //_neutrino was made with W mass constraint; use MET instead
        //      objects.push_back(m_met);
        objects.push_back(m_lepton);
        TopTopologicalVariables topo(objects);
        return topo.M();
}

double LJetsTopoVars::Jet1Jet2_M() {
  if(m_jets.size()>=2) {
    vector<TMBLorentzVector> objects;
    objects.push_back(m_jets.at(0));
    objects.push_back(m_jets.at(1));
    TopTopologicalVariables topo(objects);
    return topo.M();
  } else return -1;
}

double LJetsTopoVars::Jet1Jet2_Pt() {
  if(m_jets.size()>=2) {        
    vector<TMBLorentzVector> objects;
    objects.push_back(m_jets.at(0));
    objects.push_back(m_jets.at(1));
    TopTopologicalVariables topo(objects);
    return topo.Pt();
  } else return -1;
}

double LJetsTopoVars::Jet1Jet2_DeltaR() {
  if(m_jets.size()>=2) {  
    return m_jets.at(0).DeltaR(m_jets.at(1));
  } else return -1;
}

double LJetsTopoVars::Jet1Jet2_DeltaPhi() {
  if(m_jets.size()>=2) {  
    return TMath::Abs(m_jets.at(0).DeltaPhi(m_jets.at(1)));
  } else return -1;
}


double LJetsTopoVars::Jet1Jet2W_M() {
  if(m_jets.size()>=2) {
    vector<TMBLorentzVector> objects;
    objects.push_back(_neutrino);  //_neutrino was made with W mass constraint; use MET instead
    //  objects.push_back(m_met);
    objects.push_back(m_lepton);
    objects.push_back(m_jets.at(0));
    objects.push_back(m_jets.at(1));
    TopTopologicalVariables topo(objects);
    return topo.M();
  } else return -1;
}

double LJetsTopoVars::Jet1Jet2W_Pt() {
  if(m_jets.size()>=2) {        
    vector<TMBLorentzVector> objects;
    objects.push_back(_neutrino);  //_neutrino was made with W mass constraint; use MET instead
    //  objects.push_back(m_met);
    objects.push_back(m_lepton);
    objects.push_back(m_jets.at(0));
    objects.push_back(m_jets.at(1));
    TopTopologicalVariables topo(objects);
    return topo.Pt();
  } else return -1;
}

double LJetsTopoVars::DphiJMET() {
        return kinem::delta_phi(m_met.Phi(), m_jets.at(0).Phi());
}

double LJetsTopoVars::LeptonJet_DeltaR() {
  if(m_jets.size()>=2) {  
    return m_lepton.DeltaR(m_jets.at(0))< m_lepton.DeltaR(m_jets.at(1)) ? m_lepton.DeltaR(m_jets.at(0)) : m_lepton.DeltaR(m_jets.at(1));
  } else {
    return m_lepton.DeltaR(m_jets.at(0));
  }
}

double LJetsTopoVars::Muon_DeltaR() {
        //is this already stored in the muon somewhere?
        double DeltaR = 1e99;
        for (unsigned int i=0; i<m_jets.size(); i++) DeltaR = min(DeltaR, m_lepton.DeltaR(m_jets.at(i)));
        return DeltaR;
}


/***** removed temoprarily
double LJetsTopoVars::Muon_etHaloScaled() {
        if (!m_isMuon) throw runtime_error("LJetsTopoVars: Muon_etHaloScaled: event is not mu+jets");
        //TMBMuon* muon = dynamic_cast<TMBMuon*>(&m_lepton);
        return _muon->etHalo()/_muon->Pt();
}

double LJetsTopoVars::Muon_etTrkConeScaled() {
        if (!m_isMuon) throw runtime_error("LJetsTopoVars: Muon_etHaloScaled: event is not mu+jets");
        //TMBMuon* muon = dynamic_cast<TMBMuon*>(&m_lepton);
        return _muon->etTrkCone5()/_muon->Pt();
}

double LJetsTopoVars::Electron_iso() {
        if (m_isMuon) throw runtime_error("LJetsTopoVars: Electron_iso: event is not e+jets");
        return _electron->iso();
}

double LJetsTopoVars::Electron_lhood() {
        if (m_isMuon) throw runtime_error("LJetsTopoVars: Electron_iso: event is not e+jets");
        return _electron->Lhood8();
}
*****/

//
//_____________________________________________________________________
void LJetsTopoVars::calcHt(){
  
  //ht[ 0] = Ht 
  //ht[ 1] = Htp 
  //ht[ 2] = Htpp
  //ht[ 3] = Ht2
  //ht[ 4] = Ht2p
  //ht[ 5] = Ht2pp 
  //ht[ 6] = Ht3
  //ht[ 7] = Ht3p 
  //ht[ 8] = Ht3pp 
  //ht[ 9] = centrality
  //ht[10] = NJW;
  //ht[11] = eta_max
  //ht[12] = MdijetMin
  //ht[13] = Mtjets (definition from Jean-Roche)
  //ht[14] = sqrtsT (=Tobi's M_{T} in Note) from Tobi
  //ht[15] = MtAurelio
  //ht[16] = pZoverHT
  //ht[17] = Mevent
  //ht[18] = M123inv
  //ht[19] = Eta2Sum (Eta^2 sum)
  //ht[20] = mWrec
  //ht[21] = H = sum(jetE)

  //reset
  for(unsigned int i=0;i<_ht.size();i++) _ht[i]=0.;


  double h        = 0.;
  double hz       = 0.;
  double hx       = 0.;
  double hy       = 0.;
  double hzSigned = 0.;
  _ht[12]         =-1.;
  double mtjets   = 0.;
  TMBLorentzVector Mevent;
  int nJet = m_jets.size();
  for(int i=0;i<nJet;i++){
    hz += TMath::Abs(m_jets[i].Pz());
    hx += m_jets[i].Px();
    hy += m_jets[i].Py();
    h  += m_jets[i].E();
    _ht[0] += m_jets[i].Pt();
    Mevent += m_jets[i];
    hzSigned += m_jets[i].Pz();
    if(i>0) _ht[3] += m_jets[i].Pt();
    if(i>1) _ht[6] += m_jets[i].Pt();
    if(TMath::Abs(m_jets[i].Eta())>_ht[11] && i<4){
      _ht[11] = TMath::Abs(m_jets[i].Eta());
    }
    for(int j=i+1; j<nJet; j++){
      double mDijet = (m_jets[i]+m_jets[j]).Mag();
      if(_ht[12]<0. || mDijet<_ht[12]){ _ht[12]=mDijet; }
    }
    mtjets += 
      TMath::Power(m_jets[i].E(),2)  -
      TMath::Power(m_jets[i].Px(),2) -
      TMath::Power(m_jets[i].Py(),2);

    _ht[19] += m_jets[i].Eta()*m_jets[i].Eta();
  }
  _ht[21] = h; 

  // the "M_T"s
  if(mtjets > 0.){ _ht[13]=TMath::Sqrt(mtjets); } 
  _ht[14] = _ht[0]*_ht[0] - hx*hx - hy*hy;
  if(_ht[14]>0.) _ht[14] = TMath::Sqrt(_ht[14]);
  _ht[15] = h*h - hzSigned*hzSigned;
  if(_ht[15]>0.) _ht[15] = TMath::Sqrt(_ht[15]);

  if(_ht[0]>0.) _ht[16] = hzSigned/_ht[0];

  double hzNoLep = hz;
  hz += TMath::Abs(m_lepton.Pz());
  hz += TMath::Abs(_neutrino.Pz());
  if(hz!=0.){
    _ht[1]=_ht[0]/hz;
    _ht[4]=_ht[3]/hz;
    _ht[7]=_ht[6]/hz;
  }
  if(hzNoLep!=0.){
    _ht[2]=_ht[0]/hzNoLep;
    _ht[5]=_ht[3]/hzNoLep;
    _ht[8]=_ht[6]/hzNoLep;
  }
  if(h>0.){
    _ht[9] = _ht[0]/h;
  }

  //
  // NJW
  //
  double NJW=0;
  for(Int_t ijet=0; ijet<nJet-1; ijet++){
    double emin=55.;
    double emax=55.;
    if(m_jets[ijet  ].Pt() < 55.){emax=m_jets[ijet  ].Pt();}
    if(m_jets[ijet+1].Pt() < 55.){emin=m_jets[ijet+1].Pt();}
    NJW += 0.5*(emax*emax-emin*emin)*(ijet+1);
  }
  double elo=15.;
  if(m_jets[nJet-1].Pt() > elo){elo=m_jets[nJet-1].Pt();}
  NJW += 0.5*(elo*elo-(15.*15.))*(nJet);
  NJW /= ((55*55)-100.)/2.0;
  _ht[10] = NJW;
  

  // total event invariant mass
  Mevent += m_lepton;
  Mevent += _neutrino;
  _ht[17] = Mevent.Mag();

  // sum of dijet invariant masses for three highest jets
  // and mWrec
  if(nJet>2){
    double min=1e10;
    for(int i=0;i<2;i++){
      for(int j=i+1; j<3; j++){
        double m = (m_jets[i]+m_jets[j]).Mag();
        _ht[18] += m;
        double diff = TMath::Abs(80.4-m);
        if(diff<min){
          min = diff;
          _ht[20] = m;
        }
      }
    }
  }

  _htOK = true;
}

//
//_____________________________________________________________________
void LJetsTopoVars::calcEvtTopo(){

  //evtTopo[0] = sphericity
  //evtTopo[1] = aplanarity
  //evtTopo[2] = aplanarity including muon

  int nJet = m_jets.size();

  // calculate tensor
  //
  double psum = 0.;
  for(int k=0;k<nJet;k++){
    psum += m_jets[k].Vect().Mag32();
  }
  
  TMatrixDSym M(3);
  for(int i=0;i<3;i++){
    for(int j=i;j<3;j++){
      M(i,j)=0.;
      for(int k=0;k<nJet;k++){
        M(i,j) += m_jets[k](i) * m_jets[k](j);
      }
      M(i,j)/=psum;
      if(i!=j){M(j,i) = M(i,j);} 
    }
  }
  
  //
  // get eigenvalues
  //
  TMatrixDSymEigen eigenMatrix(M);
  const TVectorD *eigen = &eigenMatrix.GetEigenValues();

  eigenval.ResizeTo(eigen->GetNrows());
  eigenval = *eigen;
  
  //NGO fix eigenvalues to zero if too small
  //otherwise ev might be marginally below zero!
  for(int i=0;i<3;i++){
    if(fabs(eigenval[i])<1e-10) eigenval[i]=0.;
  }

  _evtTopo[0] = (3./2.) * (eigenval[1]+eigenval[2]);
  _evtTopo[1] = (3./2.) *              eigenval[2];


  //
  // some consistency checks
  //
  if( eigenval[0]<eigenval[1] || 
      eigenval[0]<eigenval[2] ||
      eigenval[1]<eigenval[2]){
    cout << "ERROR: Eigenvals not ordered!" << endl;
    std::exit(1);
  }
  if(_evtTopo[0]<0. || _evtTopo[1]<0.){
    cout << "ERROR: SPHERICITY: " << _evtTopo[0] << endl;
    cout << "ERROR: APLANARITY: " << _evtTopo[1] << endl;
    M.Print();
  }


  //
  //
  // include muon in calculation
  //
  // ------------------------------------------------------
  std::vector<TMBLorentzVector> jetMu(m_jets);
  jetMu.push_back(m_lepton);
  nJet = jetMu.size();

  // calculate tensor
  psum = 0.;
  for(int k=0;k<nJet;k++){
    psum += jetMu[k].Vect().Mag32();
  }
  
  for(int i=0;i<3;i++){
    for(int j=i;j<3;j++){
      M(i,j)=0.;
      for(int k=0;k<nJet;k++){
        M(i,j) += jetMu[k](i) * jetMu[k](j);
      }
      M(i,j)/=psum;
      if(i!=j){M(j,i) = M(i,j);} 
    }
  }
  
  // get eigenvalues
  TMatrixDSymEigen       eigenMatrix_01(M);
  TVectorD eigenval_01 = eigenMatrix_01.GetEigenValues();
    for(int i=0;i<3;i++){
    if(fabs(eigenval_01[i])<1e-10) eigenval_01[i]=0.;
  }
  _evtTopo[2] = (3./2.) * eigenval_01[2];


  _evtTopoOK = true;
}


//
//_____________________________________________________________________
void LJetsTopoVars::calcKt()
{

  //kt[0] = Ktminp
  //kt[1] = Ktminpreduced
  //kt[2] = dRmin(jet,jet);

  int nJet = m_jets.size();

  double dRmin = 9999.;
  double eTmin = 9999.;
  for(int i=0;i<nJet-1;i++){
    for(int j=i+1;j<nJet;j++){
      double dR = m_jets[i].DeltaR(m_jets[j]);
      if(dR<dRmin){
        dRmin = dR;
        eTmin = std::min(m_jets[i].Pt(),m_jets[j].Pt());
      }
    }
  }
  if(dRmin>100.) {dRmin=0.;}
  
  _kt[0] = dRmin*eTmin/(m_lepton.Pt()+_neutrino.Pt());
  _kt[1] = dRmin*eTmin;
  _kt[2] = dRmin;

  _ktOK = true;
}


//
//_____________________________________________________________________
void LJetsTopoVars::calcMt()
{
  
  //mt[0] = dPhi(muon,MET)
  //mt[1] = mT

  //
  // NGO: which met should be used in calculating the mt??
  // the raw or the one form the neutrino pz calculation, which 
  // might be scaled???
  //
  double met = TMath::Sqrt(TMath::Power(_neutrino.Px(),2.)+
                           TMath::Power(_neutrino.Py(),2));
  
    
  _mt[0] = TMath::Abs(m_lepton.DeltaPhi(_neutrino));
  _mt[1] = TMath::Sqrt(2*m_lepton.Pt()*met*(1.-TMath::Cos(_mt[0])));
  
  _mtOK = true;
}
Revision:	1.2
Committed:	Tue Feb 24 23:33:22 2009 UTC (16 years, 2 months ago) by kukartse
Content type:	text/plain
Branch:	MAIN
CVS Tags:	gak011410, gak010310, ejterm2010_25nov2009, V00-02-01, V00-02-00, gak112409, CMSSW_22X_branch_base, segala101609, V00-01-15, V00-01-14, V00-01-13, V00-01-12, V00-01-11, V00-01-10, gak031009, gak030509
Branch point for:	CMSSW_22X_branch
Changes since 1.1:	+21 -11 lines
Log Message:	* empty log message *
#	User	Rev	Content
1	kukartse	1.1	#include "LJMet/MultivariateAnalysis/interface/LJetsTopoVars.h"
2			#include "LJMet/MultivariateAnalysis/interface/TopTopologicalVariables.h"
3			#include "LJMet/MultivariateAnalysis/interface/AnglesUtil.h"
4
5			#include "TMatrixDSymEigen.h"
6
7			#include <stdexcept>
8
9			#include "DataFormats/PatCandidates/interface/Jet.h"
10			#include "DataFormats/PatCandidates/interface/Electron.h"
11			#include "DataFormats/PatCandidates/interface/Muon.h"
12			#include "DataFormats/PatCandidates/interface/MET.h"
13
14			using namespace std;
15			//using namespace cafe;
16			using namespace top_cafe;
17
18
19			namespace
20			{
21			bool moreThan(const pat::Jet& a, const pat::Jet& b)
22			{
23			return a.pt() > b.pt();
24			}
25			}
26
27
28	kukartse	1.2	int LJetsTopoVars::setEvent(const edm::Event& event, double min_dr_jet_lepton)
29	kukartse	1.1	{
30			using namespace edm;
31
32	kukartse	1.2	int removed_jets = 0;
33
34	kukartse	1.1	Handle< vector< pat::Jet > > jets;
35			Handle< vector< pat::MET > > met;
36			Handle< vector< pat::Electron > > electrons;
37			Handle< vector< pat::Muon > > muons;
38
39			event . getByLabel( m_jetBranch, jets );
40			event . getByLabel( m_metBranch, met );
41			if (m_isMuon) event . getByLabel( m_leptonBranch, muons );
42			else event . getByLabel( m_leptonBranch, electrons );
43
44			//cout << endl << "=====> LJetsTopoVars::setEvent(): number of jets = " << jets -> size() << endl;
45
46			m_jets.clear();
47
48			eigenval.ResizeTo(3);
49			eigenval.Zero();
50
51
52			if (met->size()>0){
53			m_met = TMBLorentzVector(met->begin()->pt(),met->begin()->eta(),met->begin()->phi(),met->begin()->energy(),TMBLorentzVector::kPtEtaPhiE);
54			}
55			else{
56			cout << "LJetsTopoVars::setEvent(): no MET in this event!" << endl;
57			}
58
59			if(m_isMuon){
60			if (muons->size()>0){
61			m_lepton = TMBLorentzVector(muons->begin()->pt(),muons->begin()->eta(),muons->begin()->phi(),muons->begin()->energy(),TMBLorentzVector::kPtEtaPhiE);
62			//_muon = &(event.getCollection<TMBMuon>(m_leptonBranch.c_str())[0]);
63			}
64			}
65			else{
66			if (electrons->size()>0){
67			m_lepton = TMBLorentzVector(electrons->begin()->pt(),electrons->begin()->eta(),electrons->begin()->phi(),electrons->begin()->energy(),TMBLorentzVector::kPtEtaPhiE);
68			//_electron = &(event.getCollection<TMBEMCluster>(m_leptonBranch.c_str())[0]);
69			}
70			}
71
72
73	kukartse	1.2	// loop over first four jets (that are not too close to the lepton!!!)
74			//vector<pat::Jet>::const_iterator jet;
75			//int nMax = jets->size() > 4 ? 4 : jets->size();
76			//jet = jets->begin();
77			//for (int i=0; i<nMax; ++i){
78			for (vector<pat::Jet>::const_iterator jet=jets->begin(); (jet!=jets->end()) && (m_jets.size()!=4); jet++){
79			//cout << "LJetsTopoVars::setEvent(): jet pt() = " << jet -> pt() << endl;
80			TMBLorentzVector _j(jet->pt(),jet->eta(),jet->phi(),jet->energy(),TMBLorentzVector::kPtEtaPhiE);
81			if (m_lepton.DeltaR(_j) > min_dr_jet_lepton){
82			m_jets.push_back(_j);
83			//jet++;
84			}
85			else{
86			removed_jets++;
87			}
88			}
89
90	kukartse	1.1	//original method had as given:
91			// std::vector<TheJetClass*> selectedJets,
92			// TLorentzVector selectedLepton,
93			// double nu_px (MET px OK?)
94			// double nu_py (MET py OK?)
95			double nu_px = m_met.Px();
96			double nu_py = m_met.Py();
97
98			//set all OK flags to FALSE;
99			_htOK = false;
100			_evtTopoOK = false;
101			_ktOK = false;
102			_mtOK = false;
103
104			//
105			// calculate neutrino lorentz vector (from Tobi's TopSvtAnalysis)
106			//
107			double nu_pz = 0.;
108			double nu_e = sqrt(pow(nu_px,2)+pow(nu_py,2));
109
110			double Mw = 80.43; // NGO fix this!(read from one place)
111			double l_px = m_lepton.Px();
112			double l_py = m_lepton.Py();
113			double l_pz = m_lepton.Pz();
114			double l_pt = m_lepton.Pt();
115			double l_e = m_lepton.E();
116			double Mt = sqrt(pow(l_pt+nu_e ,2)-
117			pow(l_px+nu_px,2)-
118			pow(l_py+nu_py,2));
119
120			double A;
121			if (Mt<Mw) A = pow(Mw,2)/2.;
122			else { // assume Mt=Mw, rescale MET accordingly (NGO???)
123			A = pow(Mt,2)/2.;
124			double k = nu_el_pt - nu_pxl_px - nu_py*l_py;
125			k = (k == 0. ? 0.00001 : k);
126			double scf = 0.5*pow(Mw,2)/k ;
127			nu_px *= scf;
128			nu_py *= scf;
129			nu_e = sqrt(pow(nu_px,2)+pow(nu_py,2));
130			}
131
132			double B = nu_pxl_px + nu_pyl_py;
133			double C = TMath::Max(1. + pow(nu_e,2) * (pow(l_pz,2)-pow(l_e,2)) / pow(A+B,2) , 0.);
134			C = sqrt(C);
135			double S1= (-(A+B)l_pz + (A+B)l_e*C) / (pow(l_pz,2)-pow(l_e,2));
136			double S2= (-(A+B)l_pz - (A+B)l_e*C) / (pow(l_pz,2)-pow(l_e,2));
137
138			// choose solution with smallest \|l_pz\| a la Run I
139			nu_pz = fabs (S1) < fabs (S2) ? S1 : S2 ;
140
141			//NGO: NOTE: neutrino PX, PY are not necessarily metPX, metPY any more!!!
142			_neutrino.SetPxPyPzE(nu_px,nu_py,nu_pz,nu_e);
143
144			/*****
145			*****/
146	kukartse	1.2	return removed_jets;
147	kukartse	1.1	}
148
149			double LJetsTopoVars::aplanarity() const
150			{
151			vector<TMBLorentzVector> objects(m_jets);
152			objects.push_back(m_lepton);
153			TopTopologicalVariables jetsPlusLepton(objects);
154			return jetsPlusLepton.Aplanarity();
155			}
156
157			double LJetsTopoVars::centrality() const
158			{
159			TopTopologicalVariables jets(m_jets);
160			return jets.Centrality();
161			}
162
163			double LJetsTopoVars::sphericity() const
164			{
165			vector<TMBLorentzVector> objects(m_jets);
166			objects.push_back(m_lepton);
167			TopTopologicalVariables jetsPlusLepton(objects);
168			return jetsPlusLepton.Sphericity();
169			}
170
171			double LJetsTopoVars::ht() const
172			{
173			TopTopologicalVariables jets(m_jets);
174			return jets.Ht();
175			}
176
177			double LJetsTopoVars::htpluslepton() const
178			{
179			vector<TMBLorentzVector> objects(m_jets);
180			objects.push_back(m_lepton);
181			TopTopologicalVariables jetsPlusLepton(objects);
182			return jetsPlusLepton.Ht();
183			}
184
185			double LJetsTopoVars::methtpluslepton() const
186			{
187			vector<TMBLorentzVector> objects(m_jets);
188			objects.push_back(m_lepton);
189			objects.push_back(m_met);
190			TopTopologicalVariables metjetsPlusLepton(objects);
191			return metjetsPlusLepton.Ht();
192			}
193
194			double LJetsTopoVars::h() const
195			{
196			TopTopologicalVariables jets(m_jets);
197			return jets.H();
198			}
199
200
201			double LJetsTopoVars::ktMinPrime() const
202			{
203			TopTopologicalVariables jets(m_jets);
204			float ktmin = jets.KtMin();
205			float etw = m_met.Pt() + m_lepton.Pt();
206			return ktmin/etw;
207			}
208
209			double LJetsTopoVars::dphiLepMet() const
210			{
211			return kinem::delta_phi(m_met.Phi(), m_lepton.Phi());
212			}
213
214			double LJetsTopoVars::minDijetMass() const
215			{
216			TopTopologicalVariables jets(m_jets);
217			return jets.MinimumPairMass();
218			}
219
220			double LJetsTopoVars::maxJetEta() const
221			{
222			double jetEta = 0;
223			for (unsigned int i=0; i<m_jets.size(); i++) {
224			if(TMath::Abs(m_jets.at(i).Eta()) > TMath::Abs(jetEta) ) jetEta = TMath::Abs(m_jets.at(i).Eta());
225			}
226			return jetEta;
227			}
228
229
230			double LJetsTopoVars::Et3() const
231			{
232			double Et3 = 0;
233			for (unsigned int i=2; i<m_jets.size(); i++) {
234			Et3+=m_jets.at(i).Pt();
235			}
236			return Et3;
237			}
238
239			double LJetsTopoVars::minDijetDeltaR() const
240			{
241
242			int nJet = m_jets.size();
243
244			double dRmin = 9999.;
245			double eTmin = 9999.;
246			for(int i=0;i<nJet-1;i++){
247			for(int j=i+1;j<nJet;j++){
248			double dR = m_jets[i].DeltaR(m_jets[j]);
249			if(dR<dRmin){
250			dRmin = dR;
251			eTmin = std::min(m_jets[i].Pt(),m_jets[j].Pt());
252			}
253			}
254			}
255			if(dRmin>100.) {dRmin=0.;}
256
257			return dRmin;
258			}
259
260
261			double LJetsTopoVars::Hz() {
262			vector<TMBLorentzVector> objects;
263			objects.assign(m_jets.begin(), m_jets.end());
264			objects.push_back(m_lepton);
265			objects.push_back(_neutrino);
266			double pz = 0;
267			for (vector<TMBLorentzVector>::iterator obj = objects.begin(); obj!=objects.end(); ++obj) pz += abs((*obj).Pz());
268			return pz;
269			}
270
271			double LJetsTopoVars::HT2() {
272			vector<TMBLorentzVector> objects;
273			objects.assign(++m_jets.begin(), m_jets.end());
274			TopTopologicalVariables topo(objects);
275			return topo.Ht();
276			}
277
278			double LJetsTopoVars::HT2prime() {
279			return HT2()/Hz();
280			}
281
282			double LJetsTopoVars::W_MT() {
283			vector<TMBLorentzVector> objects;
284			//objects.push_back(_neutrino); //_neutrino was made with W mass constraint; use MET instead
285			objects.push_back(m_met);
286			objects.push_back(m_lepton);
287			TopTopologicalVariables topo(objects);
288			return topo.TransverseMass();
289			}
290
291			double LJetsTopoVars::W_Pt() {
292			vector<TMBLorentzVector> objects;
293			//objects.push_back(_neutrino); //_neutrino was made with W mass constraint; use MET instead
294			objects.push_back(m_met);
295			objects.push_back(m_lepton);
296			TopTopologicalVariables topo(objects);
297			return topo.Pt();
298			}
299
300			double LJetsTopoVars::W_M() {
301			vector<TMBLorentzVector> objects;
302			objects.push_back(_neutrino); //_neutrino was made with W mass constraint; use MET instead
303			// objects.push_back(m_met);
304			objects.push_back(m_lepton);
305			TopTopologicalVariables topo(objects);
306			return topo.M();
307			}
308
309			double LJetsTopoVars::Jet1Jet2_M() {
310			if(m_jets.size()>=2) {
311			vector<TMBLorentzVector> objects;
312			objects.push_back(m_jets.at(0));
313			objects.push_back(m_jets.at(1));
314			TopTopologicalVariables topo(objects);
315			return topo.M();
316			} else return -1;
317			}
318
319			double LJetsTopoVars::Jet1Jet2_Pt() {
320			if(m_jets.size()>=2) {
321			vector<TMBLorentzVector> objects;
322			objects.push_back(m_jets.at(0));
323			objects.push_back(m_jets.at(1));
324			TopTopologicalVariables topo(objects);
325			return topo.Pt();
326			} else return -1;
327			}
328
329			double LJetsTopoVars::Jet1Jet2_DeltaR() {
330			if(m_jets.size()>=2) {
331			return m_jets.at(0).DeltaR(m_jets.at(1));
332			} else return -1;
333			}
334
335			double LJetsTopoVars::Jet1Jet2_DeltaPhi() {
336			if(m_jets.size()>=2) {
337			return TMath::Abs(m_jets.at(0).DeltaPhi(m_jets.at(1)));
338			} else return -1;
339			}
340
341
342
343			double LJetsTopoVars::Jet1Jet2W_M() {
344			if(m_jets.size()>=2) {
345			vector<TMBLorentzVector> objects;
346			objects.push_back(_neutrino); //_neutrino was made with W mass constraint; use MET instead
347			// objects.push_back(m_met);
348			objects.push_back(m_lepton);
349			objects.push_back(m_jets.at(0));
350			objects.push_back(m_jets.at(1));
351			TopTopologicalVariables topo(objects);
352			return topo.M();
353			} else return -1;
354			}
355
356			double LJetsTopoVars::Jet1Jet2W_Pt() {
357			if(m_jets.size()>=2) {
358			vector<TMBLorentzVector> objects;
359			objects.push_back(_neutrino); //_neutrino was made with W mass constraint; use MET instead
360			// objects.push_back(m_met);
361			objects.push_back(m_lepton);
362			objects.push_back(m_jets.at(0));
363			objects.push_back(m_jets.at(1));
364			TopTopologicalVariables topo(objects);
365			return topo.Pt();
366			} else return -1;
367			}
368
369			double LJetsTopoVars::DphiJMET() {
370			return kinem::delta_phi(m_met.Phi(), m_jets.at(0).Phi());
371			}
372
373			double LJetsTopoVars::LeptonJet_DeltaR() {
374			if(m_jets.size()>=2) {
375			return m_lepton.DeltaR(m_jets.at(0))< m_lepton.DeltaR(m_jets.at(1)) ? m_lepton.DeltaR(m_jets.at(0)) : m_lepton.DeltaR(m_jets.at(1));
376			} else {
377			return m_lepton.DeltaR(m_jets.at(0));
378			}
379			}
380
381			double LJetsTopoVars::Muon_DeltaR() {
382			//is this already stored in the muon somewhere?
383			double DeltaR = 1e99;
384			for (unsigned int i=0; i<m_jets.size(); i++) DeltaR = min(DeltaR, m_lepton.DeltaR(m_jets.at(i)));
385			return DeltaR;
386			}
387
388
389			/***** removed temoprarily
390			double LJetsTopoVars::Muon_etHaloScaled() {
391			if (!m_isMuon) throw runtime_error("LJetsTopoVars: Muon_etHaloScaled: event is not mu+jets");
392			//TMBMuon* muon = dynamic_cast<TMBMuon*>(&m_lepton);
393			return _muon->etHalo()/_muon->Pt();
394			}
395
396			double LJetsTopoVars::Muon_etTrkConeScaled() {
397			if (!m_isMuon) throw runtime_error("LJetsTopoVars: Muon_etHaloScaled: event is not mu+jets");
398			//TMBMuon* muon = dynamic_cast<TMBMuon*>(&m_lepton);
399			return _muon->etTrkCone5()/_muon->Pt();
400			}
401
402			double LJetsTopoVars::Electron_iso() {
403			if (m_isMuon) throw runtime_error("LJetsTopoVars: Electron_iso: event is not e+jets");
404			return _electron->iso();
405			}
406
407			double LJetsTopoVars::Electron_lhood() {
408			if (m_isMuon) throw runtime_error("LJetsTopoVars: Electron_iso: event is not e+jets");
409			return _electron->Lhood8();
410			}
411			*****/
412
413			//
414			//_____________________________________________________________________
415			void LJetsTopoVars::calcHt(){
416
417			//ht[ 0] = Ht
418			//ht[ 1] = Htp
419			//ht[ 2] = Htpp
420			//ht[ 3] = Ht2
421			//ht[ 4] = Ht2p
422			//ht[ 5] = Ht2pp
423			//ht[ 6] = Ht3
424			//ht[ 7] = Ht3p
425			//ht[ 8] = Ht3pp
426			//ht[ 9] = centrality
427			//ht[10] = NJW;
428			//ht[11] = eta_max
429			//ht[12] = MdijetMin
430			//ht[13] = Mtjets (definition from Jean-Roche)
431			//ht[14] = sqrtsT (=Tobi's M_{T} in Note) from Tobi
432			//ht[15] = MtAurelio
433			//ht[16] = pZoverHT
434			//ht[17] = Mevent
435			//ht[18] = M123inv
436			//ht[19] = Eta2Sum (Eta^2 sum)
437			//ht[20] = mWrec
438			//ht[21] = H = sum(jetE)
439
440			//reset
441			for(unsigned int i=0;i<_ht.size();i++) _ht[i]=0.;
442
443
444			double h = 0.;
445			double hz = 0.;
446			double hx = 0.;
447			double hy = 0.;
448			double hzSigned = 0.;
449			_ht[12] =-1.;
450			double mtjets = 0.;
451			TMBLorentzVector Mevent;
452			int nJet = m_jets.size();
453			for(int i=0;i<nJet;i++){
454			hz += TMath::Abs(m_jets[i].Pz());
455			hx += m_jets[i].Px();
456			hy += m_jets[i].Py();
457			h += m_jets[i].E();
458			_ht[0] += m_jets[i].Pt();
459			Mevent += m_jets[i];
460			hzSigned += m_jets[i].Pz();
461			if(i>0) _ht[3] += m_jets[i].Pt();
462			if(i>1) _ht[6] += m_jets[i].Pt();
463			if(TMath::Abs(m_jets[i].Eta())>_ht[11] && i<4){
464			_ht[11] = TMath::Abs(m_jets[i].Eta());
465			}
466			for(int j=i+1; j<nJet; j++){
467			double mDijet = (m_jets[i]+m_jets[j]).Mag();
468			if(_ht[12]<0. \|\| mDijet<_ht[12]){ _ht[12]=mDijet; }
469			}
470			mtjets +=
471			TMath::Power(m_jets[i].E(),2) -
472			TMath::Power(m_jets[i].Px(),2) -
473			TMath::Power(m_jets[i].Py(),2);
474
475			_ht[19] += m_jets[i].Eta()*m_jets[i].Eta();
476			}
477			_ht[21] = h;
478
479			// the "M_T"s
480			if(mtjets > 0.){ _ht[13]=TMath::Sqrt(mtjets); }
481			_ht[14] = _ht[0]_ht[0] - hxhx - hy*hy;
482			if(_ht[14]>0.) _ht[14] = TMath::Sqrt(_ht[14]);
483			_ht[15] = hh - hzSignedhzSigned;
484			if(_ht[15]>0.) _ht[15] = TMath::Sqrt(_ht[15]);
485
486			if(_ht[0]>0.) _ht[16] = hzSigned/_ht[0];
487
488			double hzNoLep = hz;
489			hz += TMath::Abs(m_lepton.Pz());
490			hz += TMath::Abs(_neutrino.Pz());
491			if(hz!=0.){
492			_ht[1]=_ht[0]/hz;
493			_ht[4]=_ht[3]/hz;
494			_ht[7]=_ht[6]/hz;
495			}
496			if(hzNoLep!=0.){
497			_ht[2]=_ht[0]/hzNoLep;
498			_ht[5]=_ht[3]/hzNoLep;
499			_ht[8]=_ht[6]/hzNoLep;
500			}
501			if(h>0.){
502			_ht[9] = _ht[0]/h;
503			}
504
505			//
506			// NJW
507			//
508			double NJW=0;
509			for(Int_t ijet=0; ijet<nJet-1; ijet++){
510			double emin=55.;
511			double emax=55.;
512			if(m_jets[ijet ].Pt() < 55.){emax=m_jets[ijet ].Pt();}
513			if(m_jets[ijet+1].Pt() < 55.){emin=m_jets[ijet+1].Pt();}
514			NJW += 0.5(emaxemax-eminemin)(ijet+1);
515			}
516			double elo=15.;
517			if(m_jets[nJet-1].Pt() > elo){elo=m_jets[nJet-1].Pt();}
518			NJW += 0.5(eloelo-(15.15.))(nJet);
519			NJW /= ((55*55)-100.)/2.0;
520			_ht[10] = NJW;
521
522
523			// total event invariant mass
524			Mevent += m_lepton;
525			Mevent += _neutrino;
526			_ht[17] = Mevent.Mag();
527
528			// sum of dijet invariant masses for three highest jets
529			// and mWrec
530			if(nJet>2){
531			double min=1e10;
532			for(int i=0;i<2;i++){
533			for(int j=i+1; j<3; j++){
534			double m = (m_jets[i]+m_jets[j]).Mag();
535			_ht[18] += m;
536			double diff = TMath::Abs(80.4-m);
537			if(diff<min){
538			min = diff;
539			_ht[20] = m;
540			}
541			}
542			}
543			}
544
545			_htOK = true;
546			}
547
548			//
549			//_____________________________________________________________________
550			void LJetsTopoVars::calcEvtTopo(){
551
552			//evtTopo[0] = sphericity
553			//evtTopo[1] = aplanarity
554			//evtTopo[2] = aplanarity including muon
555
556			int nJet = m_jets.size();
557
558			// calculate tensor
559			//
560			double psum = 0.;
561			for(int k=0;k<nJet;k++){
562			psum += m_jets[k].Vect().Mag32();
563			}
564
565			TMatrixDSym M(3);
566			for(int i=0;i<3;i++){
567			for(int j=i;j<3;j++){
568			M(i,j)=0.;
569			for(int k=0;k<nJet;k++){
570			M(i,j) += m_jets[k](i) * m_jets[k](j);
571			}
572			M(i,j)/=psum;
573			if(i!=j){M(j,i) = M(i,j);}
574			}
575			}
576
577			//
578			// get eigenvalues
579			//
580			TMatrixDSymEigen eigenMatrix(M);
581			const TVectorD *eigen = &eigenMatrix.GetEigenValues();
582
583			eigenval.ResizeTo(eigen->GetNrows());
584			eigenval = *eigen;
585
586			//NGO fix eigenvalues to zero if too small
587			//otherwise ev might be marginally below zero!
588			for(int i=0;i<3;i++){
589			if(fabs(eigenval[i])<1e-10) eigenval[i]=0.;
590			}
591
592			_evtTopo[0] = (3./2.) * (eigenval[1]+eigenval[2]);
593			_evtTopo[1] = (3./2.) * eigenval[2];
594
595
596			//
597			// some consistency checks
598			//
599			if( eigenval[0]<eigenval[1] \|\|
600			eigenval[0]<eigenval[2] \|\|
601			eigenval[1]<eigenval[2]){
602			cout << "ERROR: Eigenvals not ordered!" << endl;
603			std::exit(1);
604			}
605			if(_evtTopo[0]<0. \|\| _evtTopo[1]<0.){
606			cout << "ERROR: SPHERICITY: " << _evtTopo[0] << endl;
607			cout << "ERROR: APLANARITY: " << _evtTopo[1] << endl;
608			M.Print();
609			}
610
611
612			//
613			//
614			// include muon in calculation
615			//
616			// ------------------------------------------------------
617			std::vector<TMBLorentzVector> jetMu(m_jets);
618			jetMu.push_back(m_lepton);
619			nJet = jetMu.size();
620
621			// calculate tensor
622			psum = 0.;
623			for(int k=0;k<nJet;k++){
624			psum += jetMu[k].Vect().Mag32();
625			}
626
627			for(int i=0;i<3;i++){
628			for(int j=i;j<3;j++){
629			M(i,j)=0.;
630			for(int k=0;k<nJet;k++){
631			M(i,j) += jetMu[k](i) * jetMu[k](j);
632			}
633			M(i,j)/=psum;
634			if(i!=j){M(j,i) = M(i,j);}
635			}
636			}
637
638			// get eigenvalues
639			TMatrixDSymEigen eigenMatrix_01(M);
640			TVectorD eigenval_01 = eigenMatrix_01.GetEigenValues();
641			for(int i=0;i<3;i++){
642			if(fabs(eigenval_01[i])<1e-10) eigenval_01[i]=0.;
643			}
644			_evtTopo[2] = (3./2.) * eigenval_01[2];
645
646
647			_evtTopoOK = true;
648			}
649
650
651			//
652			//_____________________________________________________________________
653			void LJetsTopoVars::calcKt()
654			{
655
656			//kt[0] = Ktminp
657			//kt[1] = Ktminpreduced
658			//kt[2] = dRmin(jet,jet);
659
660			int nJet = m_jets.size();
661
662			double dRmin = 9999.;
663			double eTmin = 9999.;
664			for(int i=0;i<nJet-1;i++){
665			for(int j=i+1;j<nJet;j++){
666			double dR = m_jets[i].DeltaR(m_jets[j]);
667			if(dR<dRmin){
668			dRmin = dR;
669			eTmin = std::min(m_jets[i].Pt(),m_jets[j].Pt());
670			}
671			}
672			}
673			if(dRmin>100.) {dRmin=0.;}
674
675			_kt[0] = dRmin*eTmin/(m_lepton.Pt()+_neutrino.Pt());
676			_kt[1] = dRmin*eTmin;
677			_kt[2] = dRmin;
678
679			_ktOK = true;
680			}
681
682
683			//
684			//_____________________________________________________________________
685			void LJetsTopoVars::calcMt()
686			{
687
688			//mt[0] = dPhi(muon,MET)
689			//mt[1] = mT
690
691			//
692			// NGO: which met should be used in calculating the mt??
693			// the raw or the one form the neutrino pz calculation, which
694			// might be scaled???
695			//
696			double met = TMath::Sqrt(TMath::Power(_neutrino.Px(),2.)+
697			TMath::Power(_neutrino.Py(),2));
698
699
700			_mt[0] = TMath::Abs(m_lepton.DeltaPhi(_neutrino));
701			_mt[1] = TMath::Sqrt(2m_lepton.Pt()met*(1.-TMath::Cos(_mt[0])));
702
703			_mtOK = true;
704			}