HbbAnalysis/src/Objects.cc

#include "UserCode/HbbAnalysis/interface/Objects.hh"

namespace HbbAnalysis {

  double DeltaR(const TLorentzVector & v1, const TLorentzVector & v2){
    double dEta = v1.Eta() - v2.Eta();
    double dPhi = fabs(v1.Phi() - v2.Phi());
    if (dPhi > TMath::Pi()) dPhi = (2.0*TMath::Pi() - dPhi);
    return sqrt(dEta*dEta+dPhi*dPhi);
  }

  double DeltaPhi(const double phi1, const double phi2)
  {
    double dPhi = fabs(phi1 - phi2);
    if (dPhi > TMath::Pi()) dPhi = (2.0*TMath::Pi() - dPhi);
    //double dPhi = phi1 - phi2;
    
    return dPhi;
  }
  
  bool SameSign(double charge1, double charge2)
  {
    double product = charge1 * charge2;
    return ((product > 0.5) && (product < 1.5));
  }

  bool OppSign(double charge1, double charge2)
  {
    double product = charge1 * charge2;
    return ((product < -0.5) && (product > -1.5));
  }
  /*Fix - base vars no longer exists

  TLorentzVector FourMomentum(const BaseVars & v, const double scale) 
  {
    double lpx = v.pT*cos(v.phi);
    double lpy = v.pT*sin(v.phi);
    //double lp = v.pT/sin(2*atan(exp(-v.eta)));
    //double lpz = sqrt(lp*lp - v.pT*v.pT);
    double lpz = v.pT*sinh(v.eta);
    double lE = v.pT*cosh(v.eta);//v.E

    return TLorentzVector(lpx/scale,lpy/scale,lpz/scale,lE/scale);

  }*/

  double TransverseMass(//const BaseVars & leg1, 
                        //const BaseVars & leg2, 
                        const TLorentzVector & leg1, 
                        const TLorentzVector & leg2, 
                        const double mEx, 
                        const double mEy)
  {
    double px = leg1.Pt()*cos(leg1.Phi()) + leg2.Pt()*cos(leg2.Phi()) + mEx;
    double py = leg1.Pt()*sin(leg1.Phi()) + leg2.Pt()*sin(leg2.Phi()) + mEy;
    double et = leg1.Pt() + leg2.Pt() + TMath::Sqrt(mEx*mEx + mEy*mEy);
    double mt2 = et*et - (px*px + py*py);
    if ( mt2 < 0 ) {
      //std::cout << " --- WARNING : mt2 = " << mt2 << " is negative... Set to 0.";
      return 0.;
    }
    return sqrt(mt2);
  }

  double TransverseMass(//const BaseVars & leg1, 
                        const TLorentzVector & leg1, 
                        const double mEx, 
                        const double mEy)
  {
    double px = leg1.Pt()*cos(leg1.Phi()) + mEx;
    double py = leg1.Pt()*sin(leg1.Phi()) + mEy;
    double et = leg1.Pt() + TMath::Sqrt(mEx*mEx + mEy*mEy);
    double mt = et*et - (px*px + py*py);
    if ( mt < 0 ) {
      //std::cout << " --- WARNING : mt = " << mt << " is negative... Set to 0.";
      return 0.;
    }
    return sqrt(mt);
  }

  TLorentzVector FourMomentumCDFmethod(//const BaseVars & leg1, 
                                       //const BaseVars & leg2, 
                                       const TLorentzVector & leg1, 
                                       const TLorentzVector & leg2, 
                                       double mEx, 
                                       double mEy)
  {
    double lpx = leg1.Pt()*cos(leg1.Phi()) + leg2.Pt()*cos(leg2.Phi()) + mEx;
    double lpy = leg1.Pt()*sin(leg1.Phi()) + leg2.Pt()*sin(leg2.Phi()) + mEy;
    double lpz = leg1.Pt()*sinh(leg1.Eta()) + leg2.Pt()*sinh(leg2.Eta());
    double le =  leg1.Pt()*cosh(leg1.Eta()) + leg2.Pt()*cosh(leg2.Eta()) + TMath::Sqrt(mEx*mEx + mEy*mEy);
    return TLorentzVector(lpx, lpy, lpz, le);
  }

  TLorentzVector FourMomentumCollinearApprox(//const BaseVars & leg1,
                                             //const BaseVars & leg2,
                                             const TLorentzVector & leg1,
                                             const TLorentzVector & leg2,
                                             double mEx, 
                                             double mEy)
  {
    double px1 = leg1.Pt()*cos(leg1.Phi());
    double px2 = leg2.Pt()*cos(leg2.Phi());
    double py1 = leg1.Pt()*sin(leg1.Phi());
    double py2 = leg2.Pt()*sin(leg2.Phi());

    double x1_numerator = px1*py2 - px2*py1;
    double x1_denominator = py2*(px1 + mEx) - px2*(py1 + mEy);
    double x1 = ( x1_denominator != 0. ) ? x1_numerator/x1_denominator : -1.;

    double x2_numerator = x1_numerator;
    double x2_denominator = px1*(py2 + mEy) - py1*(px2 + mEx);
    double x2 = ( x2_denominator != 0. ) ? x2_numerator/x2_denominator : -1.;

    if ( (x1 > 0. && x1 < 1.) &&
         (x2 > 0. && x2 < 1.) ) {
      TLorentzVector p4 = leg1*x1 + leg2*x2;
      return p4;
    } else {
      return TLorentzVector(0,0,0,0);
    }
  }

  /*
  double EtaDetector(const BaseVars & v1){
    double pDet[3];
    pDet[0] = v1.pT*cos(v1.phi) + v1.vx;
    pDet[1] = v1.pT*sin(v1.phi) + v1.vy;

    double theta = 2*atan(exp(-v1.eta));
    if (pDet[1]<0) theta = TMath::Pi()+theta;

    if (tan(theta)!=0) pDet[2] = v1.pT/tan(theta) + v1.vz;
    else return -10;

    double pTDet = sqrt(pDet[0]*pDet[0] + pDet[1]*pDet[1]);
    double pDetNorm = sqrt(pDet[0]*pDet[0] + pDet[1]*pDet[1] + pDet[2]*pDet[2]);
    double thetaDet = 0;
    double cosThetaDet = 0;
    if (pDetNorm!=0) cosThetaDet = pDet[2]/pDetNorm;
    else return -10;
    if (pDet[2]!=0) thetaDet = atan(pTDet/pDet[2]);
    else return -10;
    if (cosThetaDet<0) thetaDet += TMath::Pi();
    
    return -log(tan(thetaDet/2.));
  }

  double EtaDetector(const GenVars & v1){
    double pDet[3];
    pDet[0] = v1.pT*cos(v1.phi) + v1.vx;
    pDet[1] = v1.pT*sin(v1.phi) + v1.vy;

    double theta = 2*atan(exp(-v1.eta));
    if (pDet[1]<0) theta = TMath::Pi()+theta;

    if (tan(theta)!=0) pDet[2] = v1.pT/tan(theta) + v1.vz;
    else return -10;

    double pTDet = sqrt(pDet[0]*pDet[0] + pDet[1]*pDet[1]);
    double pDetNorm = sqrt(pDet[0]*pDet[0] + pDet[1]*pDet[1] + pDet[2]*pDet[2]);
    double thetaDet = 0;
    double cosThetaDet = 0;
    if (pDetNorm!=0) cosThetaDet = pDet[2]/pDetNorm;
    else return -10;
    if (pDet[2]!=0) thetaDet = atan(pTDet/pDet[2]);
    else return -10;
    if (cosThetaDet<0) thetaDet += TMath::Pi();
    
    return -log(tan(thetaDet/2.));
  }*/


}//namespace

Revision:	1.8
Committed:	Tue Oct 25 13:26:52 2011 UTC (13 years, 6 months ago) by agilbert
Content type:	text/plain
Branch:	MAIN
Changes since 1.7:	+11 -12 lines
Log Message:	Significant code re-write. Compiles under 4_2_4 but may not work as expected. Files marked as broken may need to be fixed in the future.
#	User	Rev	Content
1	amagnan	1.1	#include "UserCode/HbbAnalysis/interface/Objects.hh"
2
3			namespace HbbAnalysis {
4	amagnan	1.7
5			double DeltaR(const TLorentzVector & v1, const TLorentzVector & v2){
6			double dEta = v1.Eta() - v2.Eta();
7			double dPhi = fabs(v1.Phi() - v2.Phi());
8			if (dPhi > TMath::Pi()) dPhi = (2.0*TMath::Pi() - dPhi);
9			return sqrt(dEtadEta+dPhidPhi);
10			}
11
12	amagnan	1.1	double DeltaPhi(const double phi1, const double phi2)
13			{
14	amagnan	1.7	double dPhi = fabs(phi1 - phi2);
15			if (dPhi > TMath::Pi()) dPhi = (2.0*TMath::Pi() - dPhi);
16			//double dPhi = phi1 - phi2;
17	amagnan	1.6
18	amagnan	1.1	return dPhi;
19			}
20	agilbert	1.8
21			bool SameSign(double charge1, double charge2)
22	amagnan	1.1	{
23	agilbert	1.8	double product = charge1 * charge2;
24			return ((product > 0.5) && (product < 1.5));
25	amagnan	1.1	}
26
27	agilbert	1.8	bool OppSign(double charge1, double charge2)
28	amagnan	1.1	{
29	agilbert	1.8	double product = charge1 * charge2;
30			return ((product < -0.5) && (product > -1.5));
31	amagnan	1.1	}
32	agilbert	1.8	/*Fix - base vars no longer exists
33	amagnan	1.1
34			TLorentzVector FourMomentum(const BaseVars & v, const double scale)
35			{
36			double lpx = v.pT*cos(v.phi);
37			double lpy = v.pT*sin(v.phi);
38			//double lp = v.pT/sin(2*atan(exp(-v.eta)));
39			//double lpz = sqrt(lplp - v.pTv.pT);
40			double lpz = v.pT*sinh(v.eta);
41			double lE = v.pT*cosh(v.eta);//v.E
42
43			return TLorentzVector(lpx/scale,lpy/scale,lpz/scale,lE/scale);
44
45	agilbert	1.8	}*/
46	amagnan	1.1
47	amagnan	1.7	double TransverseMass(//const BaseVars & leg1,
48			//const BaseVars & leg2,
49			const TLorentzVector & leg1,
50			const TLorentzVector & leg2,
51	amagnan	1.1	const double mEx,
52			const double mEy)
53			{
54	amagnan	1.7	double px = leg1.Pt()cos(leg1.Phi()) + leg2.Pt()cos(leg2.Phi()) + mEx;
55			double py = leg1.Pt()sin(leg1.Phi()) + leg2.Pt()sin(leg2.Phi()) + mEy;
56			double et = leg1.Pt() + leg2.Pt() + TMath::Sqrt(mExmEx + mEymEy);
57	amagnan	1.1	double mt2 = etet - (pxpx + py*py);
58			if ( mt2 < 0 ) {
59	amagnan	1.5	//std::cout << " --- WARNING : mt2 = " << mt2 << " is negative... Set to 0.";
60	amagnan	1.1	return 0.;
61			}
62			return sqrt(mt2);
63			}
64
65	amagnan	1.7	double TransverseMass(//const BaseVars & leg1,
66			const TLorentzVector & leg1,
67	amagnan	1.1	const double mEx,
68			const double mEy)
69			{
70	amagnan	1.7	double px = leg1.Pt()*cos(leg1.Phi()) + mEx;
71			double py = leg1.Pt()*sin(leg1.Phi()) + mEy;
72			double et = leg1.Pt() + TMath::Sqrt(mExmEx + mEymEy);
73	amagnan	1.1	double mt = etet - (pxpx + py*py);
74			if ( mt < 0 ) {
75	amagnan	1.5	//std::cout << " --- WARNING : mt = " << mt << " is negative... Set to 0.";
76	amagnan	1.1	return 0.;
77			}
78			return sqrt(mt);
79			}
80
81	amagnan	1.7	TLorentzVector FourMomentumCDFmethod(//const BaseVars & leg1,
82			//const BaseVars & leg2,
83			const TLorentzVector & leg1,
84			const TLorentzVector & leg2,
85	amagnan	1.1	double mEx,
86			double mEy)
87			{
88	amagnan	1.7	double lpx = leg1.Pt()cos(leg1.Phi()) + leg2.Pt()cos(leg2.Phi()) + mEx;
89			double lpy = leg1.Pt()sin(leg1.Phi()) + leg2.Pt()sin(leg2.Phi()) + mEy;
90			double lpz = leg1.Pt()sinh(leg1.Eta()) + leg2.Pt()sinh(leg2.Eta());
91			double le = leg1.Pt()cosh(leg1.Eta()) + leg2.Pt()cosh(leg2.Eta()) + TMath::Sqrt(mExmEx + mEymEy);
92	amagnan	1.1	return TLorentzVector(lpx, lpy, lpz, le);
93			}
94
95	amagnan	1.7	TLorentzVector FourMomentumCollinearApprox(//const BaseVars & leg1,
96			//const BaseVars & leg2,
97			const TLorentzVector & leg1,
98			const TLorentzVector & leg2,
99	amagnan	1.1	double mEx,
100			double mEy)
101			{
102	amagnan	1.7	double px1 = leg1.Pt()*cos(leg1.Phi());
103			double px2 = leg2.Pt()*cos(leg2.Phi());
104			double py1 = leg1.Pt()*sin(leg1.Phi());
105			double py2 = leg2.Pt()*sin(leg2.Phi());
106	amagnan	1.1
107			double x1_numerator = px1py2 - px2py1;
108			double x1_denominator = py2(px1 + mEx) - px2(py1 + mEy);
109			double x1 = ( x1_denominator != 0. ) ? x1_numerator/x1_denominator : -1.;
110
111			double x2_numerator = x1_numerator;
112			double x2_denominator = px1(py2 + mEy) - py1(px2 + mEx);
113			double x2 = ( x2_denominator != 0. ) ? x2_numerator/x2_denominator : -1.;
114
115			if ( (x1 > 0. && x1 < 1.) &&
116			(x2 > 0. && x2 < 1.) ) {
117	amagnan	1.7	TLorentzVector p4 = leg1x1 + leg2x2;
118	amagnan	1.1	return p4;
119			} else {
120			return TLorentzVector(0,0,0,0);
121			}
122			}
123
124	agilbert	1.8	/*
125	amagnan	1.4	double EtaDetector(const BaseVars & v1){
126			double pDet[3];
127			pDet[0] = v1.pT*cos(v1.phi) + v1.vx;
128			pDet[1] = v1.pT*sin(v1.phi) + v1.vy;
129
130			double theta = 2*atan(exp(-v1.eta));
131			if (pDet[1]<0) theta = TMath::Pi()+theta;
132
133			if (tan(theta)!=0) pDet[2] = v1.pT/tan(theta) + v1.vz;
134			else return -10;
135
136			double pTDet = sqrt(pDet[0]pDet[0] + pDet[1]pDet[1]);
137			double pDetNorm = sqrt(pDet[0]pDet[0] + pDet[1]pDet[1] + pDet[2]*pDet[2]);
138			double thetaDet = 0;
139			double cosThetaDet = 0;
140			if (pDetNorm!=0) cosThetaDet = pDet[2]/pDetNorm;
141			else return -10;
142			if (pDet[2]!=0) thetaDet = atan(pTDet/pDet[2]);
143			else return -10;
144			if (cosThetaDet<0) thetaDet += TMath::Pi();
145
146			return -log(tan(thetaDet/2.));
147			}
148
149			double EtaDetector(const GenVars & v1){
150			double pDet[3];
151			pDet[0] = v1.pT*cos(v1.phi) + v1.vx;
152			pDet[1] = v1.pT*sin(v1.phi) + v1.vy;
153
154			double theta = 2*atan(exp(-v1.eta));
155			if (pDet[1]<0) theta = TMath::Pi()+theta;
156
157			if (tan(theta)!=0) pDet[2] = v1.pT/tan(theta) + v1.vz;
158			else return -10;
159
160			double pTDet = sqrt(pDet[0]pDet[0] + pDet[1]pDet[1]);
161			double pDetNorm = sqrt(pDet[0]pDet[0] + pDet[1]pDet[1] + pDet[2]*pDet[2]);
162			double thetaDet = 0;
163			double cosThetaDet = 0;
164			if (pDetNorm!=0) cosThetaDet = pDet[2]/pDetNorm;
165			else return -10;
166			if (pDet[2]!=0) thetaDet = atan(pTDet/pDet[2]);
167			else return -10;
168			if (cosThetaDet<0) thetaDet += TMath::Pi();
169
170			return -log(tan(thetaDet/2.));
171	agilbert	1.8	}*/
172	amagnan	1.4
173
174
175	amagnan	1.1	}//namespace
176