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 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.2
Committed:	Fri Jul 23 20:24:22 2010 UTC (14 years, 9 months ago) by ceballos
Content type:	text/plain
Branch:	MAIN
Changes since 1.1:	+40 -3 lines
Log Message:	small update
#	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			Double_t JetTools::CosineOmega(const ParticleOArr *particles){
168			if(particles->GetEntries() < 2) return -999.;
169
170			TLorentzVector v_L1(particles->At(0)->Px(),particles->At(0)->Py(),particles->At(0)->Pz(),particles->At(0)->E());
171			TLorentzVector v_L2(particles->At(1)->Px(),particles->At(1)->Py(),particles->At(1)->Pz(),particles->At(1)->E());
172
173			Double_t beta = (v_L1.P()-v_L2.P())/(v_L1.Pz()-v_L2.Pz());
174
175			TVector3 B;
176			B.SetXYZ(0.0,0.0,-1.0*beta);
177
178			v_L1.Boost(B);
179			v_L2.Boost(B);
180
181			Double_t cosomega = v_L1.Vect().Dot(v_L2.Vect())/(v_L1.P()*v_L2.P());
182
183			return cosomega;
184			}
185
186			//Transverse Higgs mass
187			Double_t JetTools::MtHiggs(const CompositeParticle dilepton, const Met met, int nsel){
188			double mtHiggs = -999.0;
189			double enell,enenn,enex,eney,mll,mnu;
190
191			if (nsel == 0){ // Use of Mt mass and mnu == mll
192			enell = TMath::Sqrt(dilepton->Pt()dilepton->Pt() + dilepton->Mt()dilepton->Mt());
193			enenn = TMath::Sqrt(met->Pt() met->Pt() + dilepton->Mt()dilepton->Mt());
194			enex = dilepton->Px() + met->Px();
195			eney = dilepton->Py() + met->Py();
196			mll = dilepton->Mass();
197			mnu = mll;
198			}
199			else if(nsel == 1){ // Use of Mt mass and mnu == 0
200			enell = TMath::Sqrt(dilepton->Pt()dilepton->Pt() + dilepton->Mt()dilepton->Mt());
201			enenn = TMath::Sqrt(met->Pt() met->Pt() + 0.00.0);
202			enex = dilepton->Px() + met->Px();
203			eney = dilepton->Py() + met->Py();
204			mll = dilepton->Mass();
205			mnu = 0.0;
206			}
207			else if(nsel == 2){ // Use of M mass and mnu == mll
208			enell = TMath::Sqrt(dilepton->Pt()dilepton->Pt() + dilepton->Mass()dilepton->Mass());
209			enenn = TMath::Sqrt(met->Pt() met->Pt() + dilepton->Mass()dilepton->Mass());
210			enex = dilepton->Px() + met->Px();
211			eney = dilepton->Py() + met->Py();
212			mll = dilepton->Mass();
213			mnu = mll;
214			}
215			else if(nsel == 3){ // Use of M mass and mnu == 0
216			enell = TMath::Sqrt(dilepton->Pt()dilepton->Pt() + dilepton->Mass()dilepton->Mass());
217			enenn = TMath::Sqrt(met->Pt() met->Pt() + 0.00.0);
218			enex = dilepton->Px() + met->Px();
219			eney = dilepton->Py() + met->Py();
220			mll = dilepton->Mass();
221			mnu = 0.0;
222			}
223			else {
224			return -999.;
225			}
226
227			mtHiggs = mllmll + mnumnu + 2.0(enellenenn - enexenex - eneyeney);
228			if(mtHiggs <= 0) mtHiggs = 0.0;
229			else mtHiggs = TMath::Sqrt(mtHiggs);
230
231			return mtHiggs;
232			}