Utils/src/JetTools.cc

#include "MitPhysics/Utils/interface/JetTools.h"
#include <algorithm>
#include <vector>

ClassImp(mithep::JetTools)

using namespace mithep;

 
JetTools::JetTools()
{
   // Constructor
}

JetTools::~JetTools() 
{
  // Destructor.
}

//Remember to remove the signal from particles before inputting into the function
Double_t JetTools::NJettiness(const ParticleOArr *particles, const JetOArr *jets, double Q, double Y){
  if(particles->GetEntries() <= 0) return 0.0;

  Double_t fval = 0.0;
  Double_t fvalpart;

  for(int i=0;i<int(particles->GetEntries());i++){
    fvalpart = (particles->At(i)->Pt()) * TMath::Exp(-TMath::Abs(particles->At(i)->Eta()-Y)); 

    for(int j=0;j<int(jets->GetEntries());j++){
      fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) * 
                 (2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-particles->At(i)->Eta()))
                - 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),particles->At(i)->Phi()))));
    }
    fval = fval + fvalpart;
  }

  fval = fval / Q;

  return fval;
}

Double_t JetTools::NJettiness(const PFCandidateOArr *pfCandidates, const JetOArr *jets, double Q, double Y){
  if(pfCandidates->GetEntries() <= 0) return 0.0;

  Double_t fval = 0.0;
  Double_t fvalpart;

  for(int i=0;i<int(pfCandidates->GetEntries());i++){
    fvalpart = (pfCandidates->At(i)->Pt()) * TMath::Exp(-TMath::Abs(pfCandidates->At(i)->Eta()-Y)); 

    for(int j=0;j<int(jets->GetEntries());j++){
      fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) * 
                 (2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-pfCandidates->At(i)->Eta()))
                - 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),pfCandidates->At(i)->Phi()))));
    }
    fval = fval + fvalpart;
  }

  fval = fval / Q;

  return fval;
}

Double_t JetTools::NJettiness(const TrackOArr *tracks, const JetOArr *jets, double Q, double Y){
  if(tracks->GetEntries() <= 0) return 0.0;

  Double_t fval = 0.0;
  Double_t fvalpart;

  for(int i=0;i<int(tracks->GetEntries());i++){
    fvalpart = (tracks->At(i)->Pt()) * TMath::Exp(-TMath::Abs(tracks->At(i)->Eta()-Y));     

    for(int j=0;j<int(jets->GetEntries());j++){
      fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) * 
                 (2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-tracks->At(i)->Eta()))
                - 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),tracks->At(i)->Phi()))));
    }
    fval = fval + fvalpart;
  }

  fval = fval / Q;
  
  return fval;
}

Double_t JetTools::NJettiness(const JetOArr *jetsS, const JetOArr *jets, double Q, double Y){
  if(jetsS->GetEntries() <= 0) return 0.0;

  Double_t fval = 0.0;
  Double_t fvalpart;

  for(int i=0;i<int(jetsS->GetEntries());i++){
    fvalpart = (jetsS->At(i)->Pt()) * TMath::Exp(-TMath::Abs(jetsS->At(i)->Eta()-Y));     

    for(int j=0;j<int(jets->GetEntries());j++){
      fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) * 
                 (2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-jetsS->At(i)->Eta()))
                - 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),jetsS->At(i)->Phi()))));
    }
    fval = fval + fvalpart;
  }

  fval = fval / Q;
  
  return fval;
}

Double_t JetTools::NJettiness(const CaloTowerOArr *calos, const JetOArr *jets, double Q, double Y){
  if(calos->GetEntries() <= 0) return 0.0;

  Double_t fval = 0.0;
  Double_t fvalpart;

  for(int i=0;i<int(calos->GetEntries());i++){
    fvalpart = (calos->At(i)->Pt()) * TMath::Exp(-TMath::Abs(calos->At(i)->Eta()-Y));     

    for(int j=0;j<int(jets->GetEntries());j++){
      fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) * 
                 (2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-calos->At(i)->Eta()))
                - 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),calos->At(i)->Phi()))));
    }
    fval = fval + fvalpart;
  }

  fval = fval / Q;
  
  return fval;
}

//M_r
Double_t JetTools::M_r(const ParticleOArr *particles){

  if(particles->GetEntries() < 2) return -999.;

  Double_t E0  = particles->At(0)->E();
  Double_t E1  = particles->At(1)->E();
  Double_t Pz0 = particles->At(0)->Pz();
  Double_t Pz1 = particles->At(1)->Pz();

  Double_t den =  TMath::Power(Pz0-Pz1, 2) - TMath::Power(E0-E1,2);
  if(den <= 0) return -100.;

  return 2.0*TMath::Sqrt(TMath::Power(E0*Pz1 - E1*Pz0, 2)/den);
}

//Beta_r
Double_t JetTools::Beta_r(const ParticleOArr *particles){

  if(particles->GetEntries() < 2) return -999.;

  Double_t E0  = particles->At(0)->E();
  Double_t E1  = particles->At(1)->E();
  Double_t Pz0 = particles->At(0)->Pz();
  Double_t Pz1 = particles->At(1)->Pz();

  return (E0-E1)/(Pz0-Pz1);
}

//M_r_t
Double_t JetTools::M_r_t(const ParticleOArr *particles, const Met *met){

  if(particles->GetEntries() < 2) return -999.;

  Double_t Pt0    = particles->At(0)->Pt();
  Double_t Pt1    = particles->At(1)->Pt();
  Double_t etmiss = met->Pt();

  Double_t Px0    = particles->At(0)->Px();
  Double_t Px1    = particles->At(1)->Px();
  Double_t metx   = met->Px();
  Double_t Py0    = particles->At(0)->Py();
  Double_t Py1    = particles->At(1)->Py();
  Double_t mety   = met->Py();

  return TMath::Sqrt(0.5*etmiss*(Pt0 + Pt1) - 0.5*(metx*(Px0 + Px1) + mety*(Py0 + Py1)));
}

//Razor
Double_t JetTools::Razor(const ParticleOArr *particles, const Met *met){
  if(particles->GetEntries() < 2) return -999.;

  Double_t mr  = M_r(particles);
  Double_t mrt = M_r_t(particles,met);
  
  if(mr != 0) return mrt/mr;
  
  return -999.;
}

//Cosine Omega
Double_t JetTools::CosineOmega(const Particle *particles0, const Particle *particles1){

  TLorentzVector v_L1(particles0->Px(),particles0->Py(),particles0->Pz(),particles0->E());
  TLorentzVector v_L2(particles1->Px(),particles1->Py(),particles1->Pz(),particles1->E());

  Double_t beta = (v_L1.P()-v_L2.P())/(v_L1.Pz()-v_L2.Pz());

  TVector3 B;
  B.SetXYZ(0.0,0.0,-1.0*beta);

  v_L1.Boost(B);
  v_L2.Boost(B);

  Double_t cosomega = v_L1.Vect().Dot(v_L2.Vect())/(v_L1.P()*v_L2.P());

  return cosomega;
}

//Transverse Higgs mass
Double_t JetTools::MtHiggs(const ParticleOArr * leptons,
                           const Met *met, double metFraction[2], int nsel){
  if(leptons->Entries() < 2) return -999.0;

  double mtHiggs = -999.0;
  double enell = 0.0;
  double enenn = 0.0;
  double enex  = 0.0;
  double eney  = 0.0;
  double mll   = 0.0;
  double mnu   = 0.0;
  CompositeParticle *dilepton = new CompositeParticle();
  dilepton->AddDaughter(leptons->At(0));
  dilepton->AddDaughter(leptons->At(1));
  
  if     (nsel == 0){ // Use of Mt mass and mnu == mll
    enell = TMath::Sqrt(dilepton->Pt()*dilepton->Pt() + dilepton->Mt()*dilepton->Mt());
    enenn = TMath::Sqrt(met->Pt() *met->Pt()  + dilepton->Mt()*dilepton->Mt());
    enex  = dilepton->Px() + met->Px();
    eney  = dilepton->Py() + met->Py();
    mll   = dilepton->Mass();
    mnu   = mll;
  }
  else if(nsel == 1){ // Use of Mt mass and mnu == 0
    enell = TMath::Sqrt(dilepton->Pt()*dilepton->Pt() + dilepton->Mt()*dilepton->Mt());
    enenn = TMath::Sqrt(met->Pt() *met->Pt()  + 0.0*0.0);
    enex  = dilepton->Px() + met->Px();
    eney  = dilepton->Py() + met->Py();
    mll   = dilepton->Mass();
    mnu   = 0.0;
  }
  else if(nsel == 2){ // Use of M mass and mnu == mll
    enell = TMath::Sqrt(dilepton->Pt()*dilepton->Pt() + dilepton->Mass()*dilepton->Mass());
    enenn = TMath::Sqrt(met->Pt() *met->Pt()  + dilepton->Mass()*dilepton->Mass());
    enex  = dilepton->Px() + met->Px();
    eney  = dilepton->Py() + met->Py();
    mll   = dilepton->Mass();
    mnu   = mll;
  }
  else if(nsel == 3){ // Use of M mass and mnu == 0
    enell = TMath::Sqrt(dilepton->Pt()*dilepton->Pt() + dilepton->Mass()*dilepton->Mass());
    enenn = TMath::Sqrt(met->Pt() *met->Pt()  + 0.0*0.0);
    enex  = dilepton->Px() + met->Px();
    eney  = dilepton->Py() + met->Py();
    mll   = dilepton->Mass();
    mnu   = 0.0;
  }
  else if(nsel == 4){ // Use of Mt mass and replacing mnu using the met optimal
    enell = TMath::Sqrt(dilepton->Pt()*dilepton->Pt() + dilepton->Mt()*dilepton->Mt());
    enenn = TMath::Sqrt(met->Pt() *met->Pt()  + 0.0*0.0);
    enex  = dilepton->Px() + met->Px();
    eney  = dilepton->Py() + met->Py();
    mll   = dilepton->Mass();
    double metAuxPx[2] = {met->Px() * metFraction[0],
                          met->Px() * (1.0 - metFraction[0])};
    double metAuxPy[2] = {met->Py() * metFraction[1],
                          met->Py() * (1.0 - metFraction[1])};
    double ene = TMath::Sqrt(metAuxPx[0]*metAuxPx[0]+metAuxPy[0]*metAuxPy[0]) +
                 TMath::Sqrt(metAuxPx[1]*metAuxPx[1]+metAuxPy[1]*metAuxPy[1]);
    double px = metAuxPx[0] + metAuxPx[1];
    double py = metAuxPy[0] + metAuxPy[1];
    mnu = TMath::Sqrt(ene*ene - px*px - py*py);
  }
  else if(nsel == 5){ // Using the optimal met value
    double metAuxPx[2] = {met->Px() * metFraction[0],
                          met->Px() * (1.0 - metFraction[0])};
    double metAuxPy[2] = {met->Py() * metFraction[1],
                          met->Py() * (1.0 - metFraction[1])};
    double ene = leptons->At(0)->Pt() + leptons->At(1)->Pt() +
                 TMath::Sqrt(metAuxPx[0]*metAuxPx[0]+metAuxPy[0]*metAuxPy[0]) +
                 TMath::Sqrt(metAuxPx[1]*metAuxPx[1]+metAuxPy[1]*metAuxPy[1]);
    double px = leptons->At(0)->Px() + leptons->At(1)->Px() +
                metAuxPx[0] + metAuxPx[1];
    double py = leptons->At(0)->Py() + leptons->At(1)->Py() +
                metAuxPy[0] + metAuxPy[1];
    mtHiggs = ene*ene - px*px - py*py;
  }
  else if(nsel == 6){ // Use the formula from hep-ph:1006.4998
    mtHiggs = 2*leptons->At(0)->Pt()*leptons->At(0)->Pt() + 2*leptons->At(1)->Pt()*leptons->At(1)->Pt() + 3 * (
      leptons->At(0)->Pt()*leptons->At(1)->Pt() + met->Pt()*(leptons->At(0)->Pt()+leptons->At(1)->Pt())
      - met->Px()*dilepton->Px() - met->Py()*dilepton->Py()
      - leptons->At(0)->Px()*leptons->At(1)->Px() - leptons->At(0)->Py()*leptons->At(1)->Py());
  }
  else if(nsel == 7){ // Use of M mass and mnu == 0
    double deltaPhiDileptonMet = MathUtils::DeltaPhi(dilepton->Phi(), 
                                                     met->Phi());
    mtHiggs = 2.0*dilepton->Pt()*met->Pt()*(1.0 - cos(deltaPhiDileptonMet));
  }

  if(nsel >= 0 && nsel <= 4){
    mtHiggs = mll*mll + mnu*mnu + 2.0*(enell*enenn - enex*enex - eney*eney);
  }

  if(mtHiggs <= 0) mtHiggs = 0.0;
  else             mtHiggs = TMath::Sqrt(mtHiggs);

  delete dilepton;

  return mtHiggs;
}

Double_t JetTools::Beta(const TrackCol *tracks, Jet *jet, const Vertex *vertex, Double_t  delta_z, Double_t delta_cone){  

  if(tracks->GetEntries() <= 0) return 1.0;

  double Pt_jets_X = 0. ;
  double Pt_jets_Y = 0. ;
  double Pt_jets_X_tot = 0. ;
  double Pt_jets_Y_tot = 0. ;

  for(int i=0;i<int(tracks->GetEntries());i++){
    if(MathUtils::DeltaR(tracks->At(i)->Mom(),jet->Mom()) < delta_cone){
      Pt_jets_X_tot += tracks->At(i)->Px();
      Pt_jets_Y_tot += tracks->At(i)->Py();  
      double pDz = TMath::Abs(tracks->At(i)->DzCorrected(*vertex));
      if(pDz < delta_z){
        Pt_jets_X += tracks->At(i)->Px();
        Pt_jets_Y += tracks->At(i)->Py();   
      }
    }
  }

  if(sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot) > 0)
    return sqrt(Pt_jets_X*Pt_jets_X + Pt_jets_Y*Pt_jets_Y) / sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot);

  return 1.0;
}


Double_t JetTools::Beta(const PFJet *jet, const Vertex *vertex, Double_t  delta_z){  
  double Pt_jets= 0. ;
  double Pt_jetsTot = 0. ;

  for(UInt_t i=0;i<jet->NPFCands();i++){
    if(jet->PFCand(i)->BestTrk()){
      Pt_jetsTot += jet->PFCand(i)->BestTrk()->Pt();
      double pDz = TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertex));
      if(pDz < delta_z){
        Pt_jets += jet->PFCand(i)->BestTrk()->Pt();
      }
    }
  }

  Double_t beta = 1.0;
  if (Pt_jetsTot > 0)
    beta = Pt_jets/Pt_jetsTot;

  return beta;
}

Double_t JetTools::Beta2(const PFJet *jet, const Vertex *vertex, Double_t  delta_z){  
  double Pt_jets= 0. ;
  double Pt_jetsTot = 0. ;

  for(UInt_t i=0;i<jet->NPFCands();i++){
    if(jet->PFCand(i)->BestTrk()){
      Pt_jetsTot += jet->PFCand(i)->BestTrk()->Pt()*jet->PFCand(i)->BestTrk()->Pt();
      double pDz = TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertex));
      if(pDz < delta_z){
        Pt_jets += jet->PFCand(i)->BestTrk()->Pt()*jet->PFCand(i)->BestTrk()->Pt();
      }
    }
  }

  Double_t beta = 1.0;
  if (Pt_jetsTot > 0)
    beta = Pt_jets/Pt_jetsTot;
  return beta;
}


Bool_t  JetTools::PassBetaVertexAssociationCut(const PFJet *jet, const Vertex *referenceVertex, const VertexCol *vertices, Double_t delta_z) {

  Bool_t passBetaCut = kTRUE;
  if(vertices->GetEntries() > 0) {
    Double_t Beta = JetTools::Beta(jet, referenceVertex, 0.2);
    Double_t Beta_other = 0.0;
    for(UInt_t nv=0; nv<vertices->GetEntries(); nv++){
      if (referenceVertex == vertices->At(nv)) continue;
      Double_t BetaAux = JetTools::Beta(jet, vertices->At(nv), 0.2);
      if(BetaAux > Beta_other) Beta_other = BetaAux;
    }
    if(Beta_other > Beta) passBetaCut = kFALSE;
  }

  return passBetaCut; 

}

Bool_t  JetTools::PassBeta2VertexAssociationCut(const PFJet *jet, const Vertex *referenceVertex, const VertexCol *vertices, Double_t delta_z) {

  Bool_t passBetaCut = kTRUE;
  if(vertices->GetEntries() > 0) {
    Double_t Beta = JetTools::Beta2(jet, referenceVertex, 0.2);
    Double_t Beta_other = 0.0;
    for(UInt_t nv=0; nv<vertices->GetEntries(); nv++){
      if (referenceVertex == vertices->At(nv)) continue;
      Double_t BetaAux = JetTools::Beta2(jet, vertices->At(nv), 0.2);
      if(BetaAux > Beta_other) Beta_other = BetaAux;
    }
    if(Beta_other > Beta) passBetaCut = kFALSE;
  }

  return passBetaCut; 

}


Int_t JetTools::MaxBetaVertexIndex(const PFJet *jet, const VertexCol *vertices, Double_t  delta_z=0.2){  
  
  Int_t vertexIndex = -1;
  double beta = -0.1;
  for (UInt_t v=0; v < vertices->GetEntries(); v++){
    Double_t betaTmp = JetTools::Beta(jet, vertices->At(v), delta_z);
    if (betaTmp > beta) {
      beta = betaTmp;
      vertexIndex = v;
    }
  }
  return vertexIndex;

}

Int_t JetTools::MaxBeta2VertexIndex(const PFJet *jet, const VertexCol *vertices, Double_t  delta_z=0.2){  
  
  Int_t vertexIndex = -1;
  double beta = -0.1;
  for (UInt_t v=0; v < vertices->GetEntries(); v++){
    Double_t betaTmp = JetTools::Beta2(jet, vertices->At(v), delta_z);
    if (betaTmp > beta) {
      beta = betaTmp;
      vertexIndex = v;
    }
  }
  return vertexIndex;

}


Int_t JetTools::JetToPVAssociation(const PFJet *jet, const VertexCol *vertices, Double_t  delta_z=0.2){  

  std::vector<float> verticesPt2(vertices->GetEntries());
  for(UInt_t i=0;i<jet->NPFCands();i++){
    if(jet->PFCand(i)->BestTrk()){
      double minDZ = delta_z;
      int trackVertexIndex = -1;
      for (UInt_t v=0; v < vertices->GetEntries(); v++){
        if (minDZ > TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertices->At(v)))) {
          minDZ = TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertices->At(v)));
          trackVertexIndex = v;
        }
      }
      if (trackVertexIndex < 0) continue;
      verticesPt2[trackVertexIndex]+= jet->PFCand(i)->BestTrk()->Pt()*jet->PFCand(i)->BestTrk()->Pt();
    }
  }

  Int_t vertexIndex = 0;
  float pt2Max = 0; 
  for (uint i=0; i < verticesPt2.size(); ++i){
    if (pt2Max < verticesPt2[i]) {
      pt2Max = verticesPt2[i];
      vertexIndex = i;
    }
  }
  return vertexIndex;
}
Revision:	1.14
Committed:	Wed Apr 20 10:06:15 2011 UTC (14 years ago) by ceballos
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_025c, Mit_025b, Mit_025a, Mit_025, Mit_025pre2, Mit_024b, Mit_025pre1, Mit_024a, Mit_024, Mit_023, Mit_022a, Mit_022, Mit_021, Mit_021pre2, Mit_021pre1
Changes since 1.13:	+4 -6 lines
Log Message:	new
#	User	Rev	Content
1	ceballos	1.1	#include "MitPhysics/Utils/interface/JetTools.h"
2	mzanetti	1.13	#include <algorithm>
3			#include <vector>
4	ceballos	1.1
5			ClassImp(mithep::JetTools)
6
7			using namespace mithep;
8	mzanetti	1.13
9	ceballos	1.1
10			JetTools::JetTools()
11			{
12			// Constructor
13			}
14
15			JetTools::~JetTools()
16			{
17			// Destructor.
18			}
19
20			//Remember to remove the signal from particles before inputting into the function
21	mzanetti	1.6	Double_t JetTools::NJettiness(const ParticleOArr particles, const JetOArr jets, double Q, double Y){
22	ceballos	1.2	if(particles->GetEntries() <= 0) return 0.0;
23
24	ceballos	1.1	Double_t fval = 0.0;
25			Double_t fvalpart;
26
27			for(int i=0;i<int(particles->GetEntries());i++){
28			fvalpart = (particles->At(i)->Pt()) * TMath::Exp(-TMath::Abs(particles->At(i)->Eta()-Y));
29
30			for(int j=0;j<int(jets->GetEntries());j++){
31			fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
32			(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-particles->At(i)->Eta()))
33			- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),particles->At(i)->Phi()))));
34			}
35			fval = fval + fvalpart;
36			}
37	ceballos	1.2
38	mzanetti	1.6	fval = fval / Q;
39	ceballos	1.2
40	ceballos	1.1	return fval;
41			}
42
43	mzanetti	1.6	Double_t JetTools::NJettiness(const PFCandidateOArr pfCandidates, const JetOArr jets, double Q, double Y){
44			if(pfCandidates->GetEntries() <= 0) return 0.0;
45
46			Double_t fval = 0.0;
47			Double_t fvalpart;
48
49			for(int i=0;i<int(pfCandidates->GetEntries());i++){
50			fvalpart = (pfCandidates->At(i)->Pt()) * TMath::Exp(-TMath::Abs(pfCandidates->At(i)->Eta()-Y));
51
52			for(int j=0;j<int(jets->GetEntries());j++){
53			fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
54			(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-pfCandidates->At(i)->Eta()))
55			- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),pfCandidates->At(i)->Phi()))));
56			}
57			fval = fval + fvalpart;
58			}
59
60			fval = fval / Q;
61
62			return fval;
63			}
64
65			Double_t JetTools::NJettiness(const TrackOArr tracks, const JetOArr jets, double Q, double Y){
66	ceballos	1.2	if(tracks->GetEntries() <= 0) return 0.0;
67
68	ceballos	1.1	Double_t fval = 0.0;
69			Double_t fvalpart;
70
71			for(int i=0;i<int(tracks->GetEntries());i++){
72			fvalpart = (tracks->At(i)->Pt()) * TMath::Exp(-TMath::Abs(tracks->At(i)->Eta()-Y));
73
74			for(int j=0;j<int(jets->GetEntries());j++){
75			fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
76			(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-tracks->At(i)->Eta()))
77			- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),tracks->At(i)->Phi()))));
78			}
79			fval = fval + fvalpart;
80			}
81	ceballos	1.2
82	mzanetti	1.6	fval = fval / Q;
83	ceballos	1.2
84	ceballos	1.1	return fval;
85			}
86
87	mzanetti	1.6	Double_t JetTools::NJettiness(const JetOArr jetsS, const JetOArr jets, double Q, double Y){
88	ceballos	1.2	if(jetsS->GetEntries() <= 0) return 0.0;
89
90	ceballos	1.1	Double_t fval = 0.0;
91			Double_t fvalpart;
92
93			for(int i=0;i<int(jetsS->GetEntries());i++){
94			fvalpart = (jetsS->At(i)->Pt()) * TMath::Exp(-TMath::Abs(jetsS->At(i)->Eta()-Y));
95
96			for(int j=0;j<int(jets->GetEntries());j++){
97			fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
98			(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-jetsS->At(i)->Eta()))
99			- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),jetsS->At(i)->Phi()))));
100			}
101			fval = fval + fvalpart;
102			}
103	ceballos	1.2
104	mzanetti	1.6	fval = fval / Q;
105	ceballos	1.2
106			return fval;
107			}
108
109	mzanetti	1.6	Double_t JetTools::NJettiness(const CaloTowerOArr calos, const JetOArr jets, double Q, double Y){
110	ceballos	1.2	if(calos->GetEntries() <= 0) return 0.0;
111
112			Double_t fval = 0.0;
113			Double_t fvalpart;
114
115			for(int i=0;i<int(calos->GetEntries());i++){
116			fvalpart = (calos->At(i)->Pt()) * TMath::Exp(-TMath::Abs(calos->At(i)->Eta()-Y));
117
118			for(int j=0;j<int(jets->GetEntries());j++){
119			fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
120			(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-calos->At(i)->Eta()))
121			- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),calos->At(i)->Phi()))));
122			}
123			fval = fval + fvalpart;
124			}
125
126	mzanetti	1.6	fval = fval / Q;
127	ceballos	1.2
128	ceballos	1.1	return fval;
129			}
130
131			//M_r
132			Double_t JetTools::M_r(const ParticleOArr *particles){
133
134			if(particles->GetEntries() < 2) return -999.;
135
136			Double_t E0 = particles->At(0)->E();
137			Double_t E1 = particles->At(1)->E();
138			Double_t Pz0 = particles->At(0)->Pz();
139			Double_t Pz1 = particles->At(1)->Pz();
140
141			Double_t den = TMath::Power(Pz0-Pz1, 2) - TMath::Power(E0-E1,2);
142			if(den <= 0) return -100.;
143
144			return 2.0TMath::Sqrt(TMath::Power(E0Pz1 - E1*Pz0, 2)/den);
145			}
146
147			//Beta_r
148			Double_t JetTools::Beta_r(const ParticleOArr *particles){
149
150			if(particles->GetEntries() < 2) return -999.;
151
152			Double_t E0 = particles->At(0)->E();
153			Double_t E1 = particles->At(1)->E();
154			Double_t Pz0 = particles->At(0)->Pz();
155			Double_t Pz1 = particles->At(1)->Pz();
156
157			return (E0-E1)/(Pz0-Pz1);
158			}
159
160			//M_r_t
161			Double_t JetTools::M_r_t(const ParticleOArr particles, const Met met){
162
163			if(particles->GetEntries() < 2) return -999.;
164
165			Double_t Pt0 = particles->At(0)->Pt();
166			Double_t Pt1 = particles->At(1)->Pt();
167			Double_t etmiss = met->Pt();
168
169			Double_t Px0 = particles->At(0)->Px();
170			Double_t Px1 = particles->At(1)->Px();
171			Double_t metx = met->Px();
172			Double_t Py0 = particles->At(0)->Py();
173			Double_t Py1 = particles->At(1)->Py();
174			Double_t mety = met->Py();
175
176			return TMath::Sqrt(0.5etmiss(Pt0 + Pt1) - 0.5(metx(Px0 + Px1) + mety*(Py0 + Py1)));
177			}
178
179			//Razor
180			Double_t JetTools::Razor(const ParticleOArr particles, const Met met){
181			if(particles->GetEntries() < 2) return -999.;
182
183			Double_t mr = M_r(particles);
184			Double_t mrt = M_r_t(particles,met);
185
186			if(mr != 0) return mrt/mr;
187
188			return -999.;
189			}
190
191			//Cosine Omega
192	ceballos	1.3	Double_t JetTools::CosineOmega(const Particle particles0, const Particle particles1){
193	ceballos	1.1
194	ceballos	1.3	TLorentzVector v_L1(particles0->Px(),particles0->Py(),particles0->Pz(),particles0->E());
195			TLorentzVector v_L2(particles1->Px(),particles1->Py(),particles1->Pz(),particles1->E());
196	ceballos	1.1
197			Double_t beta = (v_L1.P()-v_L2.P())/(v_L1.Pz()-v_L2.Pz());
198
199			TVector3 B;
200			B.SetXYZ(0.0,0.0,-1.0*beta);
201
202			v_L1.Boost(B);
203			v_L2.Boost(B);
204
205			Double_t cosomega = v_L1.Vect().Dot(v_L2.Vect())/(v_L1.P()*v_L2.P());
206
207			return cosomega;
208			}
209
210			//Transverse Higgs mass
211	ceballos	1.3	Double_t JetTools::MtHiggs(const ParticleOArr * leptons,
212	ceballos	1.4	const Met *met, double metFraction[2], int nsel){
213	ceballos	1.3	if(leptons->Entries() < 2) return -999.0;
214
215	ceballos	1.1	double mtHiggs = -999.0;
216	ceballos	1.3	double enell = 0.0;
217			double enenn = 0.0;
218			double enex = 0.0;
219			double eney = 0.0;
220			double mll = 0.0;
221			double mnu = 0.0;
222			CompositeParticle *dilepton = new CompositeParticle();
223			dilepton->AddDaughter(leptons->At(0));
224			dilepton->AddDaughter(leptons->At(1));
225	ceballos	1.1
226			if (nsel == 0){ // Use of Mt mass and mnu == mll
227			enell = TMath::Sqrt(dilepton->Pt()dilepton->Pt() + dilepton->Mt()dilepton->Mt());
228			enenn = TMath::Sqrt(met->Pt() met->Pt() + dilepton->Mt()dilepton->Mt());
229			enex = dilepton->Px() + met->Px();
230			eney = dilepton->Py() + met->Py();
231			mll = dilepton->Mass();
232			mnu = mll;
233			}
234			else if(nsel == 1){ // Use of Mt mass and mnu == 0
235			enell = TMath::Sqrt(dilepton->Pt()dilepton->Pt() + dilepton->Mt()dilepton->Mt());
236			enenn = TMath::Sqrt(met->Pt() met->Pt() + 0.00.0);
237			enex = dilepton->Px() + met->Px();
238			eney = dilepton->Py() + met->Py();
239			mll = dilepton->Mass();
240			mnu = 0.0;
241			}
242			else if(nsel == 2){ // Use of M mass and mnu == mll
243			enell = TMath::Sqrt(dilepton->Pt()dilepton->Pt() + dilepton->Mass()dilepton->Mass());
244			enenn = TMath::Sqrt(met->Pt() met->Pt() + dilepton->Mass()dilepton->Mass());
245			enex = dilepton->Px() + met->Px();
246			eney = dilepton->Py() + met->Py();
247			mll = dilepton->Mass();
248			mnu = mll;
249			}
250			else if(nsel == 3){ // Use of M mass and mnu == 0
251			enell = TMath::Sqrt(dilepton->Pt()dilepton->Pt() + dilepton->Mass()dilepton->Mass());
252			enenn = TMath::Sqrt(met->Pt() met->Pt() + 0.00.0);
253			enex = dilepton->Px() + met->Px();
254			eney = dilepton->Py() + met->Py();
255			mll = dilepton->Mass();
256			mnu = 0.0;
257			}
258	ceballos	1.4	else if(nsel == 4){ // Use of Mt mass and replacing mnu using the met optimal
259			enell = TMath::Sqrt(dilepton->Pt()dilepton->Pt() + dilepton->Mt()dilepton->Mt());
260			enenn = TMath::Sqrt(met->Pt() met->Pt() + 0.00.0);
261			enex = dilepton->Px() + met->Px();
262			eney = dilepton->Py() + met->Py();
263			mll = dilepton->Mass();
264			double metAuxPx[2] = {met->Px() * metFraction[0],
265			met->Px() * (1.0 - metFraction[0])};
266			double metAuxPy[2] = {met->Py() * metFraction[1],
267			met->Py() * (1.0 - metFraction[1])};
268			double ene = TMath::Sqrt(metAuxPx[0]metAuxPx[0]+metAuxPy[0]metAuxPy[0]) +
269			TMath::Sqrt(metAuxPx[1]metAuxPx[1]+metAuxPy[1]metAuxPy[1]);
270			double px = metAuxPx[0] + metAuxPx[1];
271			double py = metAuxPy[0] + metAuxPy[1];
272			mnu = TMath::Sqrt(eneene - pxpx - py*py);
273			}
274			else if(nsel == 5){ // Using the optimal met value
275			double metAuxPx[2] = {met->Px() * metFraction[0],
276			met->Px() * (1.0 - metFraction[0])};
277			double metAuxPy[2] = {met->Py() * metFraction[1],
278			met->Py() * (1.0 - metFraction[1])};
279			double ene = leptons->At(0)->Pt() + leptons->At(1)->Pt() +
280			TMath::Sqrt(metAuxPx[0]metAuxPx[0]+metAuxPy[0]metAuxPy[0]) +
281			TMath::Sqrt(metAuxPx[1]metAuxPx[1]+metAuxPy[1]metAuxPy[1]);
282			double px = leptons->At(0)->Px() + leptons->At(1)->Px() +
283			metAuxPx[0] + metAuxPx[1];
284			double py = leptons->At(0)->Py() + leptons->At(1)->Py() +
285			metAuxPy[0] + metAuxPy[1];
286			mtHiggs = eneene - pxpx - py*py;
287			}
288			else if(nsel == 6){ // Use the formula from hep-ph:1006.4998
289	ceballos	1.3	mtHiggs = 2leptons->At(0)->Pt()leptons->At(0)->Pt() + 2leptons->At(1)->Pt()leptons->At(1)->Pt() + 3 * (
290			leptons->At(0)->Pt()leptons->At(1)->Pt() + met->Pt()(leptons->At(0)->Pt()+leptons->At(1)->Pt())
291			- met->Px()dilepton->Px() - met->Py()dilepton->Py()
292			- leptons->At(0)->Px()leptons->At(1)->Px() - leptons->At(0)->Py()leptons->At(1)->Py());
293	ceballos	1.1	}
294	ceballos	1.14	else if(nsel == 7){ // Use of M mass and mnu == 0
295			double deltaPhiDileptonMet = MathUtils::DeltaPhi(dilepton->Phi(),
296			met->Phi());
297			mtHiggs = 2.0dilepton->Pt()met->Pt()*(1.0 - cos(deltaPhiDileptonMet));
298	ceballos	1.11	}
299	ceballos	1.1
300	ceballos	1.4	if(nsel >= 0 && nsel <= 4){
301	ceballos	1.3	mtHiggs = mllmll + mnumnu + 2.0(enellenenn - enexenex - eneyeney);
302			}
303	ceballos	1.11
304	ceballos	1.1	if(mtHiggs <= 0) mtHiggs = 0.0;
305			else mtHiggs = TMath::Sqrt(mtHiggs);
306
307	ceballos	1.3	delete dilepton;
308
309	ceballos	1.1	return mtHiggs;
310			}
311	ceballos	1.5
312	ceballos	1.9	Double_t JetTools::Beta(const TrackCol tracks, Jet jet, const Vertex *vertex, Double_t delta_z, Double_t delta_cone){
313
314			if(tracks->GetEntries() <= 0) return 1.0;
315	ceballos	1.5
316			double Pt_jets_X = 0. ;
317			double Pt_jets_Y = 0. ;
318			double Pt_jets_X_tot = 0. ;
319			double Pt_jets_Y_tot = 0. ;
320
321			for(int i=0;i<int(tracks->GetEntries());i++){
322			if(MathUtils::DeltaR(tracks->At(i)->Mom(),jet->Mom()) < delta_cone){
323			Pt_jets_X_tot += tracks->At(i)->Px();
324			Pt_jets_Y_tot += tracks->At(i)->Py();
325	ceballos	1.9	double pDz = TMath::Abs(tracks->At(i)->DzCorrected(*vertex));
326	ceballos	1.5	if(pDz < delta_z){
327			Pt_jets_X += tracks->At(i)->Px();
328			Pt_jets_Y += tracks->At(i)->Py();
329			}
330			}
331			}
332
333	ceballos	1.8	if(sqrt(Pt_jets_X_totPt_jets_X_tot + Pt_jets_Y_totPt_jets_Y_tot) > 0)
334	ceballos	1.9	return sqrt(Pt_jets_XPt_jets_X + Pt_jets_YPt_jets_Y) / sqrt(Pt_jets_X_totPt_jets_X_tot + Pt_jets_Y_totPt_jets_Y_tot);
335
336			return 1.0;
337	ceballos	1.5	}
338
339	ceballos	1.7
340	ceballos	1.9	Double_t JetTools::Beta(const PFJet jet, const Vertex vertex, Double_t delta_z){
341	mzanetti	1.13	double Pt_jets= 0. ;
342			double Pt_jetsTot = 0. ;
343	ceballos	1.7
344			for(UInt_t i=0;i<jet->NPFCands();i++){
345			if(jet->PFCand(i)->BestTrk()){
346	mzanetti	1.13	Pt_jetsTot += jet->PFCand(i)->BestTrk()->Pt();
347	ceballos	1.9	double pDz = TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertex));
348	ceballos	1.7	if(pDz < delta_z){
349	mzanetti	1.13	Pt_jets += jet->PFCand(i)->BestTrk()->Pt();
350	ceballos	1.7	}
351			}
352			}
353
354	mzanetti	1.13	Double_t beta = 1.0;
355			if (Pt_jetsTot > 0)
356			beta = Pt_jets/Pt_jetsTot;
357
358			return beta;
359			}
360
361			Double_t JetTools::Beta2(const PFJet jet, const Vertex vertex, Double_t delta_z){
362			double Pt_jets= 0. ;
363			double Pt_jetsTot = 0. ;
364
365			for(UInt_t i=0;i<jet->NPFCands();i++){
366			if(jet->PFCand(i)->BestTrk()){
367			Pt_jetsTot += jet->PFCand(i)->BestTrk()->Pt()*jet->PFCand(i)->BestTrk()->Pt();
368			double pDz = TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertex));
369			if(pDz < delta_z){
370			Pt_jets += jet->PFCand(i)->BestTrk()->Pt()*jet->PFCand(i)->BestTrk()->Pt();
371			}
372			}
373			}
374	ceballos	1.9
375	mzanetti	1.13	Double_t beta = 1.0;
376			if (Pt_jetsTot > 0)
377			beta = Pt_jets/Pt_jetsTot;
378			return beta;
379	ceballos	1.7	}
380	sixie	1.10
381
382			Bool_t JetTools::PassBetaVertexAssociationCut(const PFJet jet, const Vertex referenceVertex, const VertexCol *vertices, Double_t delta_z) {
383
384			Bool_t passBetaCut = kTRUE;
385			if(vertices->GetEntries() > 0) {
386			Double_t Beta = JetTools::Beta(jet, referenceVertex, 0.2);
387			Double_t Beta_other = 0.0;
388			for(UInt_t nv=0; nv<vertices->GetEntries(); nv++){
389			if (referenceVertex == vertices->At(nv)) continue;
390			Double_t BetaAux = JetTools::Beta(jet, vertices->At(nv), 0.2);
391			if(BetaAux > Beta_other) Beta_other = BetaAux;
392			}
393			if(Beta_other > Beta) passBetaCut = kFALSE;
394			}
395
396			return passBetaCut;
397
398			}
399	mzanetti	1.13
400			Bool_t JetTools::PassBeta2VertexAssociationCut(const PFJet jet, const Vertex referenceVertex, const VertexCol *vertices, Double_t delta_z) {
401
402			Bool_t passBetaCut = kTRUE;
403			if(vertices->GetEntries() > 0) {
404			Double_t Beta = JetTools::Beta2(jet, referenceVertex, 0.2);
405			Double_t Beta_other = 0.0;
406			for(UInt_t nv=0; nv<vertices->GetEntries(); nv++){
407			if (referenceVertex == vertices->At(nv)) continue;
408			Double_t BetaAux = JetTools::Beta2(jet, vertices->At(nv), 0.2);
409			if(BetaAux > Beta_other) Beta_other = BetaAux;
410			}
411			if(Beta_other > Beta) passBetaCut = kFALSE;
412			}
413
414			return passBetaCut;
415
416			}
417
418
419			Int_t JetTools::MaxBetaVertexIndex(const PFJet jet, const VertexCol vertices, Double_t delta_z=0.2){
420
421			Int_t vertexIndex = -1;
422			double beta = -0.1;
423			for (UInt_t v=0; v < vertices->GetEntries(); v++){
424			Double_t betaTmp = JetTools::Beta(jet, vertices->At(v), delta_z);
425			if (betaTmp > beta) {
426			beta = betaTmp;
427			vertexIndex = v;
428			}
429			}
430			return vertexIndex;
431
432			}
433
434			Int_t JetTools::MaxBeta2VertexIndex(const PFJet jet, const VertexCol vertices, Double_t delta_z=0.2){
435
436			Int_t vertexIndex = -1;
437			double beta = -0.1;
438			for (UInt_t v=0; v < vertices->GetEntries(); v++){
439			Double_t betaTmp = JetTools::Beta2(jet, vertices->At(v), delta_z);
440			if (betaTmp > beta) {
441			beta = betaTmp;
442			vertexIndex = v;
443			}
444			}
445			return vertexIndex;
446
447			}
448
449
450			Int_t JetTools::JetToPVAssociation(const PFJet jet, const VertexCol vertices, Double_t delta_z=0.2){
451
452			std::vector<float> verticesPt2(vertices->GetEntries());
453			for(UInt_t i=0;i<jet->NPFCands();i++){
454			if(jet->PFCand(i)->BestTrk()){
455			double minDZ = delta_z;
456			int trackVertexIndex = -1;
457			for (UInt_t v=0; v < vertices->GetEntries(); v++){
458			if (minDZ > TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertices->At(v)))) {
459			minDZ = TMath::Abs(jet->PFCand(i)->BestTrk()->DzCorrected(*vertices->At(v)));
460			trackVertexIndex = v;
461			}
462			}
463			if (trackVertexIndex < 0) continue;
464			verticesPt2[trackVertexIndex]+= jet->PFCand(i)->BestTrk()->Pt()*jet->PFCand(i)->BestTrk()->Pt();
465			}
466			}
467
468			Int_t vertexIndex = 0;
469			float pt2Max = 0;
470			for (uint i=0; i < verticesPt2.size(); ++i){
471			if (pt2Max < verticesPt2[i]) {
472			pt2Max = verticesPt2[i];
473			vertexIndex = i;
474			}
475			}
476			return vertexIndex;
477			}