HbbAnalysis/src/Objects.cc

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

namespace HbbAnalysis {
 
  double DeltaPhi(const double phi1, const double phi2)
  {
    
    double dPhi = phi1 - phi2;
    
    return dPhi;
  }

  double SameSign(const BaseVars & v1, const BaseVars & v2)
  {
    
    return v1.charge == v2.charge;
  }

  double OppSign(const BaseVars & v1, const BaseVars & v2)
  {
    
    return (v1.charge != v2.charge && 
            v1.charge != 0 && 
            v2.charge != 0);
  }

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