HbbAnalysis/src/Objects.cc

#include "UserCode/HbbAnalysis/interface/Objects.hh"
#include <cmath>

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