Utils/src/JetTools.cc

#include "MitPhysics/Utils/interface/JetTools.h"

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, bool UseQ, Double_t 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;
  }

  if(UseQ == kTRUE) fval = fval / particles->At(0)->Pt();

  return fval;
}

Double_t JetTools::NJettiness(const TrackOArr *tracks, const JetOArr *jets, bool UseQ, Double_t 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;
  }

  if(UseQ == kTRUE) fval = fval / tracks->At(0)->Pt();
  
  return fval;
}

Double_t JetTools::NJettiness(const JetOArr *jetsS, const JetOArr *jets, bool UseQ, Double_t 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;
  }

  if(UseQ == kTRUE) fval = fval / jetsS->At(0)->Pt();
  
  return fval;
}

Double_t JetTools::NJettiness(const CaloTowerOArr *calos, const JetOArr *jets, bool UseQ, Double_t 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;
  }

  if(UseQ == kTRUE) fval = fval / calos->At(0)->Pt();
  
  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());
  }

  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;
}

void JetTools::Alpha(Double_t AlphaVar[2], const TrackCol *tracks, Jet *jet, const VertexCol *vertices, Double_t  delta_z, Double_t delta_cone){  
  AlphaVar[0] = -1.0;
  AlphaVar[1] = -1.0;
  if(tracks->GetEntries() <= 0) return;

  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(*vertices->At(0)));
      if(pDz < delta_z){
        Pt_jets_X += tracks->At(i)->Px();
        Pt_jets_Y += tracks->At(i)->Py();   
      }
    }
  }

  if(jet->Pt() > 0)
    AlphaVar[0] = sqrt(Pt_jets_X*Pt_jets_X + Pt_jets_Y*Pt_jets_Y) / jet->Pt();
  if(Pt_jets_X_tot > 0 || Pt_jets_Y_tot > 0)
    AlphaVar[1] = 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);
}

Revision:	1.5
Committed:	Mon Jan 17 17:32:49 2011 UTC (14 years, 3 months ago) by ceballos
Content type:	text/plain
Branch:	MAIN
Changes since 1.4:	+29 -0 lines
Log Message:	adding cic and alpha variables
#	User	Rev	Content
1	ceballos	1.1	#include "MitPhysics/Utils/interface/JetTools.h"
2
3			ClassImp(mithep::JetTools)
4
5			using namespace mithep;
6
7			JetTools::JetTools()
8			{
9			// Constructor
10			}
11
12			JetTools::~JetTools()
13			{
14			// Destructor.
15			}
16
17			//Remember to remove the signal from particles before inputting into the function
18	ceballos	1.2	Double_t JetTools::NJettiness(const ParticleOArr particles, const JetOArr jets, bool UseQ, Double_t Y){
19			if(particles->GetEntries() <= 0) return 0.0;
20
21	ceballos	1.1	Double_t fval = 0.0;
22			Double_t fvalpart;
23
24			for(int i=0;i<int(particles->GetEntries());i++){
25			fvalpart = (particles->At(i)->Pt()) * TMath::Exp(-TMath::Abs(particles->At(i)->Eta()-Y));
26
27			for(int j=0;j<int(jets->GetEntries());j++){
28			fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
29			(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-particles->At(i)->Eta()))
30			- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),particles->At(i)->Phi()))));
31			}
32			fval = fval + fvalpart;
33			}
34	ceballos	1.2
35			if(UseQ == kTRUE) fval = fval / particles->At(0)->Pt();
36
37	ceballos	1.1	return fval;
38			}
39
40	ceballos	1.2	Double_t JetTools::NJettiness(const TrackOArr tracks, const JetOArr jets, bool UseQ, Double_t Y){
41			if(tracks->GetEntries() <= 0) return 0.0;
42
43	ceballos	1.1	Double_t fval = 0.0;
44			Double_t fvalpart;
45
46			for(int i=0;i<int(tracks->GetEntries());i++){
47			fvalpart = (tracks->At(i)->Pt()) * TMath::Exp(-TMath::Abs(tracks->At(i)->Eta()-Y));
48
49			for(int j=0;j<int(jets->GetEntries());j++){
50			fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
51			(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-tracks->At(i)->Eta()))
52			- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),tracks->At(i)->Phi()))));
53			}
54			fval = fval + fvalpart;
55			}
56	ceballos	1.2
57			if(UseQ == kTRUE) fval = fval / tracks->At(0)->Pt();
58
59	ceballos	1.1	return fval;
60			}
61
62	ceballos	1.2	Double_t JetTools::NJettiness(const JetOArr jetsS, const JetOArr jets, bool UseQ, Double_t Y){
63			if(jetsS->GetEntries() <= 0) return 0.0;
64
65	ceballos	1.1	Double_t fval = 0.0;
66			Double_t fvalpart;
67
68			for(int i=0;i<int(jetsS->GetEntries());i++){
69			fvalpart = (jetsS->At(i)->Pt()) * TMath::Exp(-TMath::Abs(jetsS->At(i)->Eta()-Y));
70
71			for(int j=0;j<int(jets->GetEntries());j++){
72			fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
73			(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-jetsS->At(i)->Eta()))
74			- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),jetsS->At(i)->Phi()))));
75			}
76			fval = fval + fvalpart;
77			}
78	ceballos	1.2
79			if(UseQ == kTRUE) fval = fval / jetsS->At(0)->Pt();
80
81			return fval;
82			}
83
84			Double_t JetTools::NJettiness(const CaloTowerOArr calos, const JetOArr jets, bool UseQ, Double_t Y){
85			if(calos->GetEntries() <= 0) return 0.0;
86
87			Double_t fval = 0.0;
88			Double_t fvalpart;
89
90			for(int i=0;i<int(calos->GetEntries());i++){
91			fvalpart = (calos->At(i)->Pt()) * TMath::Exp(-TMath::Abs(calos->At(i)->Eta()-Y));
92
93			for(int j=0;j<int(jets->GetEntries());j++){
94			fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
95			(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-calos->At(i)->Eta()))
96			- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),calos->At(i)->Phi()))));
97			}
98			fval = fval + fvalpart;
99			}
100
101			if(UseQ == kTRUE) fval = fval / calos->At(0)->Pt();
102
103	ceballos	1.1	return fval;
104			}
105
106			//M_r
107			Double_t JetTools::M_r(const ParticleOArr *particles){
108
109			if(particles->GetEntries() < 2) return -999.;
110
111			Double_t E0 = particles->At(0)->E();
112			Double_t E1 = particles->At(1)->E();
113			Double_t Pz0 = particles->At(0)->Pz();
114			Double_t Pz1 = particles->At(1)->Pz();
115
116			Double_t den = TMath::Power(Pz0-Pz1, 2) - TMath::Power(E0-E1,2);
117			if(den <= 0) return -100.;
118
119			return 2.0TMath::Sqrt(TMath::Power(E0Pz1 - E1*Pz0, 2)/den);
120			}
121
122			//Beta_r
123			Double_t JetTools::Beta_r(const ParticleOArr *particles){
124
125			if(particles->GetEntries() < 2) return -999.;
126
127			Double_t E0 = particles->At(0)->E();
128			Double_t E1 = particles->At(1)->E();
129			Double_t Pz0 = particles->At(0)->Pz();
130			Double_t Pz1 = particles->At(1)->Pz();
131
132			return (E0-E1)/(Pz0-Pz1);
133			}
134
135			//M_r_t
136			Double_t JetTools::M_r_t(const ParticleOArr particles, const Met met){
137
138			if(particles->GetEntries() < 2) return -999.;
139
140			Double_t Pt0 = particles->At(0)->Pt();
141			Double_t Pt1 = particles->At(1)->Pt();
142			Double_t etmiss = met->Pt();
143
144			Double_t Px0 = particles->At(0)->Px();
145			Double_t Px1 = particles->At(1)->Px();
146			Double_t metx = met->Px();
147			Double_t Py0 = particles->At(0)->Py();
148			Double_t Py1 = particles->At(1)->Py();
149			Double_t mety = met->Py();
150
151			return TMath::Sqrt(0.5etmiss(Pt0 + Pt1) - 0.5(metx(Px0 + Px1) + mety*(Py0 + Py1)));
152			}
153
154			//Razor
155			Double_t JetTools::Razor(const ParticleOArr particles, const Met met){
156			if(particles->GetEntries() < 2) return -999.;
157
158			Double_t mr = M_r(particles);
159			Double_t mrt = M_r_t(particles,met);
160
161			if(mr != 0) return mrt/mr;
162
163			return -999.;
164			}
165
166			//Cosine Omega
167	ceballos	1.3	Double_t JetTools::CosineOmega(const Particle particles0, const Particle particles1){
168	ceballos	1.1
169	ceballos	1.3	TLorentzVector v_L1(particles0->Px(),particles0->Py(),particles0->Pz(),particles0->E());
170			TLorentzVector v_L2(particles1->Px(),particles1->Py(),particles1->Pz(),particles1->E());
171	ceballos	1.1
172			Double_t beta = (v_L1.P()-v_L2.P())/(v_L1.Pz()-v_L2.Pz());
173
174			TVector3 B;
175			B.SetXYZ(0.0,0.0,-1.0*beta);
176
177			v_L1.Boost(B);
178			v_L2.Boost(B);
179
180			Double_t cosomega = v_L1.Vect().Dot(v_L2.Vect())/(v_L1.P()*v_L2.P());
181
182			return cosomega;
183			}
184
185			//Transverse Higgs mass
186	ceballos	1.3	Double_t JetTools::MtHiggs(const ParticleOArr * leptons,
187	ceballos	1.4	const Met *met, double metFraction[2], int nsel){
188	ceballos	1.3	if(leptons->Entries() < 2) return -999.0;
189
190	ceballos	1.1	double mtHiggs = -999.0;
191	ceballos	1.3	double enell = 0.0;
192			double enenn = 0.0;
193			double enex = 0.0;
194			double eney = 0.0;
195			double mll = 0.0;
196			double mnu = 0.0;
197			CompositeParticle *dilepton = new CompositeParticle();
198			dilepton->AddDaughter(leptons->At(0));
199			dilepton->AddDaughter(leptons->At(1));
200	ceballos	1.1
201			if (nsel == 0){ // Use of Mt mass and mnu == mll
202			enell = TMath::Sqrt(dilepton->Pt()dilepton->Pt() + dilepton->Mt()dilepton->Mt());
203			enenn = TMath::Sqrt(met->Pt() met->Pt() + dilepton->Mt()dilepton->Mt());
204			enex = dilepton->Px() + met->Px();
205			eney = dilepton->Py() + met->Py();
206			mll = dilepton->Mass();
207			mnu = mll;
208			}
209			else if(nsel == 1){ // Use of Mt mass and mnu == 0
210			enell = TMath::Sqrt(dilepton->Pt()dilepton->Pt() + dilepton->Mt()dilepton->Mt());
211			enenn = TMath::Sqrt(met->Pt() met->Pt() + 0.00.0);
212			enex = dilepton->Px() + met->Px();
213			eney = dilepton->Py() + met->Py();
214			mll = dilepton->Mass();
215			mnu = 0.0;
216			}
217			else if(nsel == 2){ // Use of M mass and mnu == mll
218			enell = TMath::Sqrt(dilepton->Pt()dilepton->Pt() + dilepton->Mass()dilepton->Mass());
219			enenn = TMath::Sqrt(met->Pt() met->Pt() + dilepton->Mass()dilepton->Mass());
220			enex = dilepton->Px() + met->Px();
221			eney = dilepton->Py() + met->Py();
222			mll = dilepton->Mass();
223			mnu = mll;
224			}
225			else if(nsel == 3){ // Use of M mass and mnu == 0
226			enell = TMath::Sqrt(dilepton->Pt()dilepton->Pt() + dilepton->Mass()dilepton->Mass());
227			enenn = TMath::Sqrt(met->Pt() met->Pt() + 0.00.0);
228			enex = dilepton->Px() + met->Px();
229			eney = dilepton->Py() + met->Py();
230			mll = dilepton->Mass();
231			mnu = 0.0;
232			}
233	ceballos	1.4	else if(nsel == 4){ // Use of Mt mass and replacing mnu using the met optimal
234			enell = TMath::Sqrt(dilepton->Pt()dilepton->Pt() + dilepton->Mt()dilepton->Mt());
235			enenn = TMath::Sqrt(met->Pt() met->Pt() + 0.00.0);
236			enex = dilepton->Px() + met->Px();
237			eney = dilepton->Py() + met->Py();
238			mll = dilepton->Mass();
239			double metAuxPx[2] = {met->Px() * metFraction[0],
240			met->Px() * (1.0 - metFraction[0])};
241			double metAuxPy[2] = {met->Py() * metFraction[1],
242			met->Py() * (1.0 - metFraction[1])};
243			double ene = TMath::Sqrt(metAuxPx[0]metAuxPx[0]+metAuxPy[0]metAuxPy[0]) +
244			TMath::Sqrt(metAuxPx[1]metAuxPx[1]+metAuxPy[1]metAuxPy[1]);
245			double px = metAuxPx[0] + metAuxPx[1];
246			double py = metAuxPy[0] + metAuxPy[1];
247			mnu = TMath::Sqrt(eneene - pxpx - py*py);
248			}
249			else if(nsel == 5){ // Using the optimal met value
250			double metAuxPx[2] = {met->Px() * metFraction[0],
251			met->Px() * (1.0 - metFraction[0])};
252			double metAuxPy[2] = {met->Py() * metFraction[1],
253			met->Py() * (1.0 - metFraction[1])};
254			double ene = leptons->At(0)->Pt() + leptons->At(1)->Pt() +
255			TMath::Sqrt(metAuxPx[0]metAuxPx[0]+metAuxPy[0]metAuxPy[0]) +
256			TMath::Sqrt(metAuxPx[1]metAuxPx[1]+metAuxPy[1]metAuxPy[1]);
257			double px = leptons->At(0)->Px() + leptons->At(1)->Px() +
258			metAuxPx[0] + metAuxPx[1];
259			double py = leptons->At(0)->Py() + leptons->At(1)->Py() +
260			metAuxPy[0] + metAuxPy[1];
261			mtHiggs = eneene - pxpx - py*py;
262			}
263			else if(nsel == 6){ // Use the formula from hep-ph:1006.4998
264	ceballos	1.3	mtHiggs = 2leptons->At(0)->Pt()leptons->At(0)->Pt() + 2leptons->At(1)->Pt()leptons->At(1)->Pt() + 3 * (
265			leptons->At(0)->Pt()leptons->At(1)->Pt() + met->Pt()(leptons->At(0)->Pt()+leptons->At(1)->Pt())
266			- met->Px()dilepton->Px() - met->Py()dilepton->Py()
267			- leptons->At(0)->Px()leptons->At(1)->Px() - leptons->At(0)->Py()leptons->At(1)->Py());
268	ceballos	1.1	}
269
270	ceballos	1.4	if(nsel >= 0 && nsel <= 4){
271	ceballos	1.3	mtHiggs = mllmll + mnumnu + 2.0(enellenenn - enexenex - eneyeney);
272			}
273	ceballos	1.1	if(mtHiggs <= 0) mtHiggs = 0.0;
274			else mtHiggs = TMath::Sqrt(mtHiggs);
275
276	ceballos	1.3	delete dilepton;
277
278	ceballos	1.1	return mtHiggs;
279			}
280	ceballos	1.5
281			void JetTools::Alpha(Double_t AlphaVar[2], const TrackCol tracks, Jet jet, const VertexCol *vertices, Double_t delta_z, Double_t delta_cone){
282			AlphaVar[0] = -1.0;
283			AlphaVar[1] = -1.0;
284			if(tracks->GetEntries() <= 0) return;
285
286			double Pt_jets_X = 0. ;
287			double Pt_jets_Y = 0. ;
288			double Pt_jets_X_tot = 0. ;
289			double Pt_jets_Y_tot = 0. ;
290
291			for(int i=0;i<int(tracks->GetEntries());i++){
292			if(MathUtils::DeltaR(tracks->At(i)->Mom(),jet->Mom()) < delta_cone){
293			Pt_jets_X_tot += tracks->At(i)->Px();
294			Pt_jets_Y_tot += tracks->At(i)->Py();
295			double pDz = TMath::Abs(tracks->At(i)->DzCorrected(*vertices->At(0)));
296			if(pDz < delta_z){
297			Pt_jets_X += tracks->At(i)->Px();
298			Pt_jets_Y += tracks->At(i)->Py();
299			}
300			}
301			}
302
303			if(jet->Pt() > 0)
304			AlphaVar[0] = sqrt(Pt_jets_XPt_jets_X + Pt_jets_YPt_jets_Y) / jet->Pt();
305			if(Pt_jets_X_tot > 0 \|\| Pt_jets_Y_tot > 0)
306			AlphaVar[1] = 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);
307			}
308