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, Double_t Y){
  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;
  }
  return fval;
}

Double_t JetTools::NJettiness(const TrackOArr *tracks, const JetOArr *jets, Double_t Y){
  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;
  }
  return fval;
}

Double_t JetTools::NJettiness(const JetOArr *jetsS, const JetOArr *jets, Double_t Y){
  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;
  }
  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 ParticleOArr *particles){
  if(particles->GetEntries() < 2) return -999.;

  TLorentzVector v_L1(particles->At(0)->Px(),particles->At(0)->Py(),particles->At(0)->Pz(),particles->At(0)->E());
  TLorentzVector v_L2(particles->At(1)->Px(),particles->At(1)->Py(),particles->At(1)->Pz(),particles->At(1)->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 CompositeParticle *dilepton, const Met *met, int nsel){
  double mtHiggs = -999.0;
  double enell,enenn,enex,eney,mll,mnu;
  
  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 {
    return -999.;
  }

  mtHiggs = mll*mll + mnu*mnu + 2.0*(enell*enenn - enex*enex - eney*eney);
  if(mtHiggs <= 0) mtHiggs = 0.0;
  else             mtHiggs = TMath::Sqrt(mtHiggs);

  return mtHiggs;
}
Revision:	1.1
Committed:	Sun Jul 18 21:14:47 2010 UTC (14 years, 9 months ago) by ceballos
Content type:	text/plain
Branch:	MAIN
Log Message:	adding jet tools
#	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			Double_t JetTools::NJettiness(const ParticleOArr particles, const JetOArr jets, Double_t Y){
19			Double_t fval = 0.0;
20			Double_t fvalpart;
21
22			for(int i=0;i<int(particles->GetEntries());i++){
23			fvalpart = (particles->At(i)->Pt()) * TMath::Exp(-TMath::Abs(particles->At(i)->Eta()-Y));
24
25			for(int j=0;j<int(jets->GetEntries());j++){
26			fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
27			(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-particles->At(i)->Eta()))
28			- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),particles->At(i)->Phi()))));
29			}
30			fval = fval + fvalpart;
31			}
32			return fval;
33			}
34
35			Double_t JetTools::NJettiness(const TrackOArr tracks, const JetOArr jets, Double_t Y){
36			Double_t fval = 0.0;
37			Double_t fvalpart;
38
39			for(int i=0;i<int(tracks->GetEntries());i++){
40			fvalpart = (tracks->At(i)->Pt()) * TMath::Exp(-TMath::Abs(tracks->At(i)->Eta()-Y));
41
42			for(int j=0;j<int(jets->GetEntries());j++){
43			fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
44			(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-tracks->At(i)->Eta()))
45			- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),tracks->At(i)->Phi()))));
46			}
47			fval = fval + fvalpart;
48			}
49			return fval;
50			}
51
52			Double_t JetTools::NJettiness(const JetOArr jetsS, const JetOArr jets, Double_t Y){
53			Double_t fval = 0.0;
54			Double_t fvalpart;
55
56			for(int i=0;i<int(jetsS->GetEntries());i++){
57			fvalpart = (jetsS->At(i)->Pt()) * TMath::Exp(-TMath::Abs(jetsS->At(i)->Eta()-Y));
58
59			for(int j=0;j<int(jets->GetEntries());j++){
60			fvalpart = TMath::Min(fvalpart,(jets->At(j)->Pt()) *
61			(2 * TMath::CosH(TMath::Abs(jets->At(j)->Eta()-jetsS->At(i)->Eta()))
62			- 2 * TMath::Cos(MathUtils::DeltaPhi(jets->At(j)->Phi(),jetsS->At(i)->Phi()))));
63			}
64			fval = fval + fvalpart;
65			}
66			return fval;
67			}
68
69			//M_r
70			Double_t JetTools::M_r(const ParticleOArr *particles){
71
72			if(particles->GetEntries() < 2) return -999.;
73
74			Double_t E0 = particles->At(0)->E();
75			Double_t E1 = particles->At(1)->E();
76			Double_t Pz0 = particles->At(0)->Pz();
77			Double_t Pz1 = particles->At(1)->Pz();
78
79			Double_t den = TMath::Power(Pz0-Pz1, 2) - TMath::Power(E0-E1,2);
80			if(den <= 0) return -100.;
81
82			return 2.0TMath::Sqrt(TMath::Power(E0Pz1 - E1*Pz0, 2)/den);
83			}
84
85			//Beta_r
86			Double_t JetTools::Beta_r(const ParticleOArr *particles){
87
88			if(particles->GetEntries() < 2) return -999.;
89
90			Double_t E0 = particles->At(0)->E();
91			Double_t E1 = particles->At(1)->E();
92			Double_t Pz0 = particles->At(0)->Pz();
93			Double_t Pz1 = particles->At(1)->Pz();
94
95			return (E0-E1)/(Pz0-Pz1);
96			}
97
98			//M_r_t
99			Double_t JetTools::M_r_t(const ParticleOArr particles, const Met met){
100
101			if(particles->GetEntries() < 2) return -999.;
102
103			Double_t Pt0 = particles->At(0)->Pt();
104			Double_t Pt1 = particles->At(1)->Pt();
105			Double_t etmiss = met->Pt();
106
107			Double_t Px0 = particles->At(0)->Px();
108			Double_t Px1 = particles->At(1)->Px();
109			Double_t metx = met->Px();
110			Double_t Py0 = particles->At(0)->Py();
111			Double_t Py1 = particles->At(1)->Py();
112			Double_t mety = met->Py();
113
114			return TMath::Sqrt(0.5etmiss(Pt0 + Pt1) - 0.5(metx(Px0 + Px1) + mety*(Py0 + Py1)));
115			}
116
117			//Razor
118			Double_t JetTools::Razor(const ParticleOArr particles, const Met met){
119			if(particles->GetEntries() < 2) return -999.;
120
121			Double_t mr = M_r(particles);
122			Double_t mrt = M_r_t(particles,met);
123
124			if(mr != 0) return mrt/mr;
125
126			return -999.;
127			}
128
129			//Cosine Omega
130			Double_t JetTools::CosineOmega(const ParticleOArr *particles){
131			if(particles->GetEntries() < 2) return -999.;
132
133			TLorentzVector v_L1(particles->At(0)->Px(),particles->At(0)->Py(),particles->At(0)->Pz(),particles->At(0)->E());
134			TLorentzVector v_L2(particles->At(1)->Px(),particles->At(1)->Py(),particles->At(1)->Pz(),particles->At(1)->E());
135
136			Double_t beta = (v_L1.P()-v_L2.P())/(v_L1.Pz()-v_L2.Pz());
137
138			TVector3 B;
139			B.SetXYZ(0.0,0.0,-1.0*beta);
140
141			v_L1.Boost(B);
142			v_L2.Boost(B);
143
144			Double_t cosomega = v_L1.Vect().Dot(v_L2.Vect())/(v_L1.P()*v_L2.P());
145
146			return cosomega;
147			}
148
149			//Transverse Higgs mass
150			Double_t JetTools::MtHiggs(const CompositeParticle dilepton, const Met met, int nsel){
151			double mtHiggs = -999.0;
152			double enell,enenn,enex,eney,mll,mnu;
153
154			if (nsel == 0){ // Use of Mt mass and mnu == mll
155			enell = TMath::Sqrt(dilepton->Pt()dilepton->Pt() + dilepton->Mt()dilepton->Mt());
156			enenn = TMath::Sqrt(met->Pt() met->Pt() + dilepton->Mt()dilepton->Mt());
157			enex = dilepton->Px() + met->Px();
158			eney = dilepton->Py() + met->Py();
159			mll = dilepton->Mass();
160			mnu = mll;
161			}
162			else if(nsel == 1){ // Use of Mt mass and mnu == 0
163			enell = TMath::Sqrt(dilepton->Pt()dilepton->Pt() + dilepton->Mt()dilepton->Mt());
164			enenn = TMath::Sqrt(met->Pt() met->Pt() + 0.00.0);
165			enex = dilepton->Px() + met->Px();
166			eney = dilepton->Py() + met->Py();
167			mll = dilepton->Mass();
168			mnu = 0.0;
169			}
170			else if(nsel == 2){ // Use of M mass and mnu == mll
171			enell = TMath::Sqrt(dilepton->Pt()dilepton->Pt() + dilepton->Mass()dilepton->Mass());
172			enenn = TMath::Sqrt(met->Pt() met->Pt() + dilepton->Mass()dilepton->Mass());
173			enex = dilepton->Px() + met->Px();
174			eney = dilepton->Py() + met->Py();
175			mll = dilepton->Mass();
176			mnu = mll;
177			}
178			else if(nsel == 3){ // Use of M mass and mnu == 0
179			enell = TMath::Sqrt(dilepton->Pt()dilepton->Pt() + dilepton->Mass()dilepton->Mass());
180			enenn = TMath::Sqrt(met->Pt() met->Pt() + 0.00.0);
181			enex = dilepton->Px() + met->Px();
182			eney = dilepton->Py() + met->Py();
183			mll = dilepton->Mass();
184			mnu = 0.0;
185			}
186			else {
187			return -999.;
188			}
189
190			mtHiggs = mllmll + mnumnu + 2.0(enellenenn - enexenex - eneyeney);
191			if(mtHiggs <= 0) mtHiggs = 0.0;
192			else mtHiggs = TMath::Sqrt(mtHiggs);
193
194			return mtHiggs;
195			}