SFrameTools/src/Utils.cxx

#include "../include/Utils.h"


//variable HEP Tagger from Rebekka

bool variableHepTopTag(TopJet topjet, double ptJetMin, double massWindowLower, double massWindowUpper, double cutCondition2, double cutCondition3){
 
  double mjet;
  double ptjet;
  int nsubjets;

  double topmass=172.3;
  double wmass=80.4;

  nsubjets=topjet.numberOfDaughters();

  LorentzVector allsubjets(0,0,0,0);
  
  for(int j=0; j<topjet.numberOfDaughters(); ++j){
    allsubjets += topjet.subjets()[j].v4();
  }
  if(!allsubjets.isTimelike()){
    mjet=0;
    return false;
  }

  mjet = allsubjets.M();
  ptjet= allsubjets.Pt();

  double m12, m13, m23;

  //The subjetcs have to be three
  if(nsubjets==3) {

    std::vector<Particle> subjets = topjet.subjets();
    sort(subjets.begin(), subjets.end(), HigherPt());

    m12 = 0;
    if( (subjets[0].v4()+subjets[1].v4()).isTimelike()) 
      m12=(subjets[0].v4()+subjets[1].v4()).M();
    m13 = 0;
    if( (subjets[0].v4()+subjets[2].v4()).isTimelike() ) 
      m13=(subjets[0].v4()+subjets[2].v4()).M();
    m23 = 0;
    if( (subjets[1].v4()+subjets[2].v4()).isTimelike()  )
      m23 = (subjets[1].v4()+subjets[2].v4()).M();
   
  }
  else{
    return false;
  }

  double rmin=massWindowLower*wmass/topmass;
  double rmax=massWindowUpper*wmass/topmass;

  int keep=0;

  //Conditions on the subjects: at least one has to be true
  //1 condition
  if(atan(m13/m12)>0.2 && atan(m13/m12)<1.3 && m23/mjet>rmin && m23/mjet<rmax) keep=1;
  
  //2 condition
  double cond2left=pow(rmin,2)*(1+pow((m13/m12),2));
  double cond2cent=1-pow(m23/mjet,2);
  double cond2right=pow(rmax,2)*(1+pow(m13/m12,2));

  if(cond2left<cond2cent && cond2cent<cond2right && m23/mjet>cutCondition2) keep=1;

  //3 condition
  double cond3left=pow(rmin,2)*(1+pow((m12/m13),2));
  double cond3cent=1-pow(m23/mjet,2);
  double cond3right=pow(rmax,2)*(1+pow(m12/m13,2));

  if(cond3left<cond3cent && cond3cent<cond3right && m23/mjet>cutCondition3) keep=1;

  //Final requirement: at least one of the three subjets conditions and total pt
  if(keep==1 && ptjet>ptJetMin){
    return true;
  }
  else{
    return false;
  }

}


//HEP Tagger from Ivan

bool HepTopTag(TopJet topjet){
 
  double mjet;
  double ptjet;
  int nsubjets;

  double topmass=172.3;
  double wmass=80.4;

  nsubjets=topjet.numberOfDaughters();

  LorentzVector allsubjets(0,0,0,0);
  
  for(int j=0; j<topjet.numberOfDaughters(); ++j){
    allsubjets += topjet.subjets()[j].v4();
  }
  if(!allsubjets.isTimelike()){
    mjet=0;
    return false;
  }

  mjet = allsubjets.M();
  ptjet= allsubjets.Pt();

  double m12, m13, m23;

  //The subjetcs have to be three
  if(nsubjets==3) {

    std::vector<Particle> subjets = topjet.subjets();
    sort(subjets.begin(), subjets.end(), HigherPt());

    m12 = 0;
    if( (subjets[0].v4()+subjets[1].v4()).isTimelike()) 
      m12=(subjets[0].v4()+subjets[1].v4()).M();
    m13 = 0;
    if( (subjets[0].v4()+subjets[2].v4()).isTimelike() ) 
      m13=(subjets[0].v4()+subjets[2].v4()).M();
    m23 = 0;
    if( (subjets[1].v4()+subjets[2].v4()).isTimelike()  )
      m23 = (subjets[1].v4()+subjets[2].v4()).M();
   
  }
  else{
    return false;
  }

  double rmin=0.85*wmass/topmass;
  double rmax=1.15*wmass/topmass;

  int keep=0;

  //Conditions on the subjects: at least one has to be true
  //1 condition
  if(atan(m13/m12)>0.2 && atan(m13/m12)<1.3 && m23/mjet>rmin && m23/mjet<rmax) keep=1;
  
  //2 condition
  double cond2left=pow(rmin,2)*(1+pow((m13/m12),2));
  double cond2cent=1-pow(m23/mjet,2);
  double cond2right=pow(rmax,2)*(1+pow(m13/m12,2));

  if(cond2left<cond2cent && cond2cent<cond2right && m23/mjet>0.35) keep=1;

  //3 condition
  double cond3left=pow(rmin,2)*(1+pow((m12/m13),2));
  double cond3cent=1-pow(m23/mjet,2);
  double cond3right=pow(rmax,2)*(1+pow(m12/m13,2));

  if(cond3left<cond3cent && cond3cent<cond3right && m23/mjet>0.35) keep=1;

  //Final requirement: at least one of the three subjets conditions and total pt
  if(keep==1 && ptjet>200){
    return true;
  }
  else{
    return false;
  }

}


//default values (mminLower=50., mjetLower=140, mjetUpper=250.) defined in Utils.h
 bool variableTopTag(TopJet topjet, double &mjet, int &nsubjets, double &mmin, double mminLower, double mjetLower, double mjetUpper){
 
  nsubjets=topjet.numberOfDaughters();
 
  LorentzVector allsubjets(0,0,0,0);
 
  for(int j=0; j<topjet.numberOfDaughters(); ++j){
  allsubjets += topjet.subjets()[j].v4();
  }
  if(!allsubjets.isTimelike()){
  mjet=0;
  mmin=0;
 //  mminLower=50;
//   mjetLower=140;
//   mjetUpper=250;
  return false;
  }
 
  mjet = allsubjets.M();
 
  if(nsubjets>=3) {
 
  std::vector<Particle> subjets = topjet.subjets();
  sort(subjets.begin(), subjets.end(), HigherPt());
 
  double m01 = 0;
  if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
  m01=(subjets[0].v4()+subjets[1].v4()).M();
  double m02 = 0;
  if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
  m02=(subjets[0].v4()+subjets[2].v4()).M();
  double m12 = 0;
  if( (subjets[1].v4()+subjets[2].v4()).isTimelike() )
  m12 = (subjets[1].v4()+subjets[2].v4()).M();
 
  //minimum pairwise mass
  mmin = std::min(m01,std::min(m02,m12));
  }
 
  //at least 3 sub-jets
  if(nsubjets<3) return false;
  //minimum pairwise mass > 50 GeV/c^2
  if(mmin<mminLower) return false;
  //jet mass between 140 and 250 GeV/c^2
  if(mjet<mjetLower || mjet>mjetUpper) return false;
 
  return true;
 }


bool TopTag(TopJet topjet,  double &mjet, int &nsubjets, double &mmin){
 
  nsubjets=topjet.numberOfDaughters();

  LorentzVector allsubjets(0,0,0,0);
  
  for(int j=0; j<topjet.numberOfDaughters(); ++j){
    allsubjets += topjet.subjets()[j].v4();
  }
  if(!allsubjets.isTimelike()){
    mjet=0;
    mmin=0;
    return false;
  }

  mjet = allsubjets.M();
  
  if(nsubjets>=3) {

    std::vector<Particle> subjets = topjet.subjets();
    sort(subjets.begin(), subjets.end(), HigherPt());

    double m01 = 0;
    if( (subjets[0].v4()+subjets[1].v4()).isTimelike()) 
      m01=(subjets[0].v4()+subjets[1].v4()).M();
    double m02 = 0;
    if( (subjets[0].v4()+subjets[2].v4()).isTimelike() ) 
      m02=(subjets[0].v4()+subjets[2].v4()).M();
    double m12 = 0;
    if( (subjets[1].v4()+subjets[2].v4()).isTimelike()  )
      m12 = (subjets[1].v4()+subjets[2].v4()).M();
    
    //minimum pairwise mass 
    mmin = std::min(m01,std::min(m02,m12));
  }

  //at least 3 sub-jets
  if(nsubjets<3) return false;
  //minimum pairwise mass > 50 GeV/c^2
  if(mmin<50) return false;
  //jet mass between 140 and 250 GeV/c^2
  if(mjet<140 || mjet>250) return false;
  
  return true;
}

Jet* nextJet(const Particle *p, std::vector<Jet> *jets){

  double deltarmin = double_infinity();
  Jet* nextjet=0;
  for(unsigned int i=0; i<jets->size();++i){
    if(jets->at(i).deltaR(*p) < deltarmin){
      deltarmin = jets->at(i).deltaR(*p);
      nextjet = &jets->at(i);
    }
  }

  return nextjet;
}

bool WTag(TopJet prunedjet,  double& mjet, int &nsubjets, double& massdrop){

  nsubjets=prunedjet.numberOfDaughters();

  mjet = 0;
  if(prunedjet.v4().isTimelike())
    mjet = prunedjet.v4().M(); 

  //calculate mass drop for first sub-jet ordered in pt
  massdrop = 0;
  if(nsubjets>=1 && mjet>0){

    std::vector< Particle > subjets = prunedjet.subjets();
    sort(subjets.begin(), subjets.end(), HigherPt());

    double m1 = 0;
    if(subjets[0].v4().isTimelike())
      m1 = subjets[0].v4().M();

    massdrop = m1/mjet;
  }
  
  //at least 2 sub-jets
  if(nsubjets<2) return false;
  //60 GeV < pruned jet mass < 100 GeV
  if(mjet <= 60 || mjet >= 100) return false;
  //mass drop < 0.4
  if(massdrop>=0.4) return false;

  return true;

}

double pTrel(const Particle *p, std::vector<Jet> *jets){

  double ptrel=0;

  Jet* nextjet =  nextJet(p,jets);

  TVector3 p3(p->v4().Px(),p->v4().Py(),p->v4().Pz());
  TVector3 jet3(nextjet->v4().Px(),nextjet->v4().Py(),nextjet->v4().Pz());

  if(p3.Mag()!=0 && jet3.Mag()!=0){
    double sin_alpha = (p3.Cross(jet3)).Mag()/p3.Mag()/jet3.Mag();
    ptrel = p3.Mag()*sin_alpha;
  }
  else{
    std::cout << "something strange happend in the ptrel calculation: either lepton or jet momentum is 0" <<std::endl;
  }

  return ptrel;
}

double deltaRmin(const Particle *p, std::vector<Jet> *jets){
  return nextJet(p,jets)->deltaR(*p);
}

TVector3 toVector(LorentzVector v4){

  TVector3 v3(0,0,0);
  v3.SetX(v4.X());
  v3.SetY(v4.Y());
  v3.SetZ(v4.Z());
  return v3;
}

TVector3 toVector(LorentzVectorXYZE v4){

  TVector3 v3(0,0,0);
  v3.SetX(v4.X());
  v3.SetY(v4.Y());
  v3.SetZ(v4.Z());
  return v3;
}

LorentzVectorXYZE toXYZ(LorentzVector v4){

  LorentzVectorXYZE v4_new(0,0,0,0);
  v4_new.SetPx(v4.Px());
  v4_new.SetPy(v4.Py());
  v4_new.SetPz(v4.Pz());
  v4_new.SetE(v4.E());
  return v4_new;
}

LorentzVector toPtEtaPhi(LorentzVectorXYZE v4){

  LorentzVector v4_new(0,0,0,0);
  v4_new.SetPt(v4.Pt());
  v4_new.SetEta(v4.Eta());
  v4_new.SetPhi(v4.Phi());
  v4_new.SetE(v4.E());
  return v4_new;
}

double deltaR(LorentzVector v1, LorentzVector v2){

  Particle p1;
  p1.set_v4(v1);
  Particle p2;
  p2.set_v4(v2);
  return p1.deltaR(p2);
}

double double_infinity(){
  return std::numeric_limits<double>::infinity() ;
}

int int_infinity(){
  return std::numeric_limits<int>::max() ;
}

int myPow(int x, unsigned int p) {
  int i = 1;
  for (unsigned int j = 1; j <= p; j++)  i *= x;
  return i;
}
Revision:	1.6
Committed:	Mon Oct 15 15:35:51 2012 UTC (12 years, 6 months ago) by hoeing
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Nov-30-2012-v2, Nov-30-2012-v1
Changes since 1.5:	+221 -0 lines
Log Message:	HEP top tagging tools added
#	User	Rev	Content
1	peiffer	1.1	#include "../include/Utils.h"
2
3
4	hoeing	1.6	//variable HEP Tagger from Rebekka
5
6			bool variableHepTopTag(TopJet topjet, double ptJetMin, double massWindowLower, double massWindowUpper, double cutCondition2, double cutCondition3){
7
8			double mjet;
9			double ptjet;
10			int nsubjets;
11
12			double topmass=172.3;
13			double wmass=80.4;
14
15			nsubjets=topjet.numberOfDaughters();
16
17			LorentzVector allsubjets(0,0,0,0);
18
19			for(int j=0; j<topjet.numberOfDaughters(); ++j){
20			allsubjets += topjet.subjets()[j].v4();
21			}
22			if(!allsubjets.isTimelike()){
23			mjet=0;
24			return false;
25			}
26
27			mjet = allsubjets.M();
28			ptjet= allsubjets.Pt();
29
30			double m12, m13, m23;
31
32			//The subjetcs have to be three
33			if(nsubjets==3) {
34
35			std::vector<Particle> subjets = topjet.subjets();
36			sort(subjets.begin(), subjets.end(), HigherPt());
37
38			m12 = 0;
39			if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
40			m12=(subjets[0].v4()+subjets[1].v4()).M();
41			m13 = 0;
42			if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
43			m13=(subjets[0].v4()+subjets[2].v4()).M();
44			m23 = 0;
45			if( (subjets[1].v4()+subjets[2].v4()).isTimelike() )
46			m23 = (subjets[1].v4()+subjets[2].v4()).M();
47
48			}
49			else{
50			return false;
51			}
52
53			double rmin=massWindowLower*wmass/topmass;
54			double rmax=massWindowUpper*wmass/topmass;
55
56			int keep=0;
57
58			//Conditions on the subjects: at least one has to be true
59			//1 condition
60			if(atan(m13/m12)>0.2 && atan(m13/m12)<1.3 && m23/mjet>rmin && m23/mjet<rmax) keep=1;
61
62			//2 condition
63			double cond2left=pow(rmin,2)*(1+pow((m13/m12),2));
64			double cond2cent=1-pow(m23/mjet,2);
65			double cond2right=pow(rmax,2)*(1+pow(m13/m12,2));
66
67			if(cond2left<cond2cent && cond2cent<cond2right && m23/mjet>cutCondition2) keep=1;
68
69			//3 condition
70			double cond3left=pow(rmin,2)*(1+pow((m12/m13),2));
71			double cond3cent=1-pow(m23/mjet,2);
72			double cond3right=pow(rmax,2)*(1+pow(m12/m13,2));
73
74			if(cond3left<cond3cent && cond3cent<cond3right && m23/mjet>cutCondition3) keep=1;
75
76			//Final requirement: at least one of the three subjets conditions and total pt
77			if(keep==1 && ptjet>ptJetMin){
78			return true;
79			}
80			else{
81			return false;
82			}
83
84			}
85
86
87
88
89
90
91			//HEP Tagger from Ivan
92
93			bool HepTopTag(TopJet topjet){
94
95			double mjet;
96			double ptjet;
97			int nsubjets;
98
99			double topmass=172.3;
100			double wmass=80.4;
101
102			nsubjets=topjet.numberOfDaughters();
103
104			LorentzVector allsubjets(0,0,0,0);
105
106			for(int j=0; j<topjet.numberOfDaughters(); ++j){
107			allsubjets += topjet.subjets()[j].v4();
108			}
109			if(!allsubjets.isTimelike()){
110			mjet=0;
111			return false;
112			}
113
114			mjet = allsubjets.M();
115			ptjet= allsubjets.Pt();
116
117			double m12, m13, m23;
118
119			//The subjetcs have to be three
120			if(nsubjets==3) {
121
122			std::vector<Particle> subjets = topjet.subjets();
123			sort(subjets.begin(), subjets.end(), HigherPt());
124
125			m12 = 0;
126			if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
127			m12=(subjets[0].v4()+subjets[1].v4()).M();
128			m13 = 0;
129			if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
130			m13=(subjets[0].v4()+subjets[2].v4()).M();
131			m23 = 0;
132			if( (subjets[1].v4()+subjets[2].v4()).isTimelike() )
133			m23 = (subjets[1].v4()+subjets[2].v4()).M();
134
135			}
136			else{
137			return false;
138			}
139
140			double rmin=0.85*wmass/topmass;
141			double rmax=1.15*wmass/topmass;
142
143			int keep=0;
144
145			//Conditions on the subjects: at least one has to be true
146			//1 condition
147			if(atan(m13/m12)>0.2 && atan(m13/m12)<1.3 && m23/mjet>rmin && m23/mjet<rmax) keep=1;
148
149			//2 condition
150			double cond2left=pow(rmin,2)*(1+pow((m13/m12),2));
151			double cond2cent=1-pow(m23/mjet,2);
152			double cond2right=pow(rmax,2)*(1+pow(m13/m12,2));
153
154			if(cond2left<cond2cent && cond2cent<cond2right && m23/mjet>0.35) keep=1;
155
156			//3 condition
157			double cond3left=pow(rmin,2)*(1+pow((m12/m13),2));
158			double cond3cent=1-pow(m23/mjet,2);
159			double cond3right=pow(rmax,2)*(1+pow(m12/m13,2));
160
161			if(cond3left<cond3cent && cond3cent<cond3right && m23/mjet>0.35) keep=1;
162
163			//Final requirement: at least one of the three subjets conditions and total pt
164			if(keep==1 && ptjet>200){
165			return true;
166			}
167			else{
168			return false;
169			}
170
171			}
172
173
174			//default values (mminLower=50., mjetLower=140, mjetUpper=250.) defined in Utils.h
175			bool variableTopTag(TopJet topjet, double &mjet, int &nsubjets, double &mmin, double mminLower, double mjetLower, double mjetUpper){
176
177			nsubjets=topjet.numberOfDaughters();
178
179			LorentzVector allsubjets(0,0,0,0);
180
181			for(int j=0; j<topjet.numberOfDaughters(); ++j){
182			allsubjets += topjet.subjets()[j].v4();
183			}
184			if(!allsubjets.isTimelike()){
185			mjet=0;
186			mmin=0;
187			// mminLower=50;
188			// mjetLower=140;
189			// mjetUpper=250;
190			return false;
191			}
192
193			mjet = allsubjets.M();
194
195			if(nsubjets>=3) {
196
197			std::vector<Particle> subjets = topjet.subjets();
198			sort(subjets.begin(), subjets.end(), HigherPt());
199
200			double m01 = 0;
201			if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
202			m01=(subjets[0].v4()+subjets[1].v4()).M();
203			double m02 = 0;
204			if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
205			m02=(subjets[0].v4()+subjets[2].v4()).M();
206			double m12 = 0;
207			if( (subjets[1].v4()+subjets[2].v4()).isTimelike() )
208			m12 = (subjets[1].v4()+subjets[2].v4()).M();
209
210			//minimum pairwise mass
211			mmin = std::min(m01,std::min(m02,m12));
212			}
213
214			//at least 3 sub-jets
215			if(nsubjets<3) return false;
216			//minimum pairwise mass > 50 GeV/c^2
217			if(mmin<mminLower) return false;
218			//jet mass between 140 and 250 GeV/c^2
219			if(mjet<mjetLower \|\| mjet>mjetUpper) return false;
220
221			return true;
222			}
223
224
225	peiffer	1.1	bool TopTag(TopJet topjet, double &mjet, int &nsubjets, double &mmin){
226
227			nsubjets=topjet.numberOfDaughters();
228
229			LorentzVector allsubjets(0,0,0,0);
230
231			for(int j=0; j<topjet.numberOfDaughters(); ++j){
232			allsubjets += topjet.subjets()[j].v4();
233			}
234			if(!allsubjets.isTimelike()){
235			mjet=0;
236			mmin=0;
237			return false;
238			}
239
240			mjet = allsubjets.M();
241
242			if(nsubjets>=3) {
243
244			std::vector<Particle> subjets = topjet.subjets();
245			sort(subjets.begin(), subjets.end(), HigherPt());
246
247			double m01 = 0;
248			if( (subjets[0].v4()+subjets[1].v4()).isTimelike())
249			m01=(subjets[0].v4()+subjets[1].v4()).M();
250			double m02 = 0;
251			if( (subjets[0].v4()+subjets[2].v4()).isTimelike() )
252			m02=(subjets[0].v4()+subjets[2].v4()).M();
253			double m12 = 0;
254			if( (subjets[1].v4()+subjets[2].v4()).isTimelike() )
255			m12 = (subjets[1].v4()+subjets[2].v4()).M();
256
257	peiffer	1.5	//minimum pairwise mass
258	peiffer	1.1	mmin = std::min(m01,std::min(m02,m12));
259			}
260
261			//at least 3 sub-jets
262			if(nsubjets<3) return false;
263			//minimum pairwise mass > 50 GeV/c^2
264			if(mmin<50) return false;
265			//jet mass between 140 and 250 GeV/c^2
266			if(mjet<140 \|\| mjet>250) return false;
267
268			return true;
269			}
270	peiffer	1.2
271			Jet* nextJet(const Particle p, std::vector<Jet> jets){
272
273			double deltarmin = double_infinity();
274			Jet* nextjet=0;
275			for(unsigned int i=0; i<jets->size();++i){
276			if(jets->at(i).deltaR(*p) < deltarmin){
277			deltarmin = jets->at(i).deltaR(*p);
278			nextjet = &jets->at(i);
279			}
280			}
281
282			return nextjet;
283			}
284
285	peiffer	1.5	bool WTag(TopJet prunedjet, double& mjet, int &nsubjets, double& massdrop){
286
287			nsubjets=prunedjet.numberOfDaughters();
288
289			mjet = 0;
290			if(prunedjet.v4().isTimelike())
291			mjet = prunedjet.v4().M();
292
293			//calculate mass drop for first sub-jet ordered in pt
294			massdrop = 0;
295			if(nsubjets>=1 && mjet>0){
296
297			std::vector< Particle > subjets = prunedjet.subjets();
298			sort(subjets.begin(), subjets.end(), HigherPt());
299
300			double m1 = 0;
301			if(subjets[0].v4().isTimelike())
302			m1 = subjets[0].v4().M();
303
304			massdrop = m1/mjet;
305			}
306
307			//at least 2 sub-jets
308			if(nsubjets<2) return false;
309			//60 GeV < pruned jet mass < 100 GeV
310			if(mjet <= 60 \|\| mjet >= 100) return false;
311			//mass drop < 0.4
312			if(massdrop>=0.4) return false;
313
314			return true;
315
316			}
317
318	peiffer	1.2	double pTrel(const Particle p, std::vector<Jet> jets){
319
320			double ptrel=0;
321
322			Jet* nextjet = nextJet(p,jets);
323
324			TVector3 p3(p->v4().Px(),p->v4().Py(),p->v4().Pz());
325			TVector3 jet3(nextjet->v4().Px(),nextjet->v4().Py(),nextjet->v4().Pz());
326
327			if(p3.Mag()!=0 && jet3.Mag()!=0){
328			double sin_alpha = (p3.Cross(jet3)).Mag()/p3.Mag()/jet3.Mag();
329			ptrel = p3.Mag()*sin_alpha;
330			}
331			else{
332			std::cout << "something strange happend in the ptrel calculation: either lepton or jet momentum is 0" <<std::endl;
333			}
334
335			return ptrel;
336			}
337
338			double deltaRmin(const Particle p, std::vector<Jet> jets){
339			return nextJet(p,jets)->deltaR(*p);
340			}
341
342	peiffer	1.4	TVector3 toVector(LorentzVector v4){
343
344			TVector3 v3(0,0,0);
345			v3.SetX(v4.X());
346			v3.SetY(v4.Y());
347			v3.SetZ(v4.Z());
348			return v3;
349			}
350
351			TVector3 toVector(LorentzVectorXYZE v4){
352
353			TVector3 v3(0,0,0);
354			v3.SetX(v4.X());
355			v3.SetY(v4.Y());
356			v3.SetZ(v4.Z());
357			return v3;
358			}
359
360			LorentzVectorXYZE toXYZ(LorentzVector v4){
361
362			LorentzVectorXYZE v4_new(0,0,0,0);
363			v4_new.SetPx(v4.Px());
364			v4_new.SetPy(v4.Py());
365			v4_new.SetPz(v4.Pz());
366			v4_new.SetE(v4.E());
367			return v4_new;
368			}
369
370			LorentzVector toPtEtaPhi(LorentzVectorXYZE v4){
371
372			LorentzVector v4_new(0,0,0,0);
373			v4_new.SetPt(v4.Pt());
374			v4_new.SetEta(v4.Eta());
375			v4_new.SetPhi(v4.Phi());
376			v4_new.SetE(v4.E());
377			return v4_new;
378			}
379
380			double deltaR(LorentzVector v1, LorentzVector v2){
381
382			Particle p1;
383			p1.set_v4(v1);
384			Particle p2;
385			p2.set_v4(v2);
386			return p1.deltaR(p2);
387			}
388	peiffer	1.2
389			double double_infinity(){
390			return std::numeric_limits<double>::infinity() ;
391			}
392
393			int int_infinity(){
394			return std::numeric_limits<int>::max() ;
395			}
396	peiffer	1.4
397			int myPow(int x, unsigned int p) {
398			int i = 1;
399			for (unsigned int j = 1; j <= p; j++) i *= x;
400			return i;
401			}