MultivariateAnalysis/interface/AnglesUtil.h

#ifndef INC_ANGLESUTIL
#define INC_ANGLESUTIL
///////////////////////////////////////////////////////////////////////////////
// File: AnglesUtil.hpp
// 
// Purpose:  Provide useful functions for calculating angles and eta
//
// Created:   4-NOV-1998  Serban Protopopescu
// History:   replaced old KinemUtil  
// Modified:  23-July-2000 Add rapidity calculation
//            14-Aug-2003 converted all functions to double precision
///////////////////////////////////////////////////////////////////////////////
// Dependencies (#includes)

#include <cmath>
#include <cstdlib>


namespace kinem{
  const double PI=2.0*acos(0.);
  const double TWOPI=2.0*PI;
  const float ETA_LIMIT=15.0;
  const float EPSILON=1.E-10;

  //  calculate phi from x, y
  inline double phi(double x, double y);
  //  calculate phi for a line defined by xy1 and xy2 (xy2-xy1)
  inline double phi(double xy1[2], double xy2[2]);
  inline double phi(float xy1[2], float xy2[2]);

  //  calculate theta from x, y, z
  inline double theta(double x, double y, double z);
  //  calculate theta for a line defined by xyz1 and xyz2 (xyz2-xyz1)
  inline double theta(double xyz1[3], double xyz2[3]);
  inline double theta(float xyz1[3], float xyz2[3]);
  //  calculate theta from eta
  inline double theta(double etap);

  //  calculate eta from x, y, z (return also theta)
  inline double eta(double x, double y, double z);
  //  calculate eta for a line defined by xyz1 and xyz2 (xyz2-xyz1)
  inline double eta(double xyz1[3], double xyz2[3]);
  inline double eta(float xyz1[3], float xyz2[3]);
  //  calculate eta from theta
  inline double eta(double th);

  // calculate rapidity from E, pz
  inline double y(double E, double pz);

  //  calculate phi1-phi2 keeping value between 0 and pi
  inline double delta_phi(double ph11, double phi2);
    // calculate phi1-phi2 keeping value between -pi and pi
  inline double signed_delta_phi(double ph11, double phi2);
  //  calculate deltaR
  inline double delta_R(double eta1, double phi1, double eta2, double phi2);

  //  calculate unit vectors given two points
  inline void uvectors(double u[3], double xyz1[3], double xyz2[3]);
  inline void uvectors(float u[3], float xyz1[3], float xyz2[3]);

  inline double tanl_from_theta(double theta);
  inline double theta_from_tanl(double tanl); 
}

inline
double kinem::phi(double x, double y)
{
  double PHI=atan2(y, x);
  return (PHI>=0)? PHI : kinem::TWOPI+PHI;
}

inline 
double kinem::phi(float xy1[2], float xy2[2]){
    double dxy1[2]={xy1[0],xy1[1]};
    double dxy2[2]={xy2[0],xy2[1]};
    return phi(dxy1,dxy2);
  }

inline
double kinem::phi(double xy1[2], double xy2[2])
{
  double x=xy2[0]-xy1[0];
  double y=xy2[1]-xy1[1];
  return phi(x, y);
}

inline
double kinem::delta_phi(double phi1, double phi2)
{
  double PHI=fabs(phi1-phi2);
  return (PHI<=PI)? PHI : kinem::TWOPI-PHI;
}

inline
double kinem::signed_delta_phi(double phi1, double phi2)
{
  double phia=phi1;
  if(phi1>PI) phia=phi1-kinem::TWOPI;
  double phib=phi2;
  if(phi2>PI) phib=phi2-kinem::TWOPI;
  double dphi=phia-phib;
  if(dphi>PI) dphi-=kinem::TWOPI;
  if(dphi<-PI) dphi+=kinem::TWOPI;
  return dphi;
}

inline double kinem::delta_R(double eta1, double phi1, double eta2, double phi2)
{
   double deta = eta1-eta2;
   double dphi = kinem::delta_phi(phi1,phi2);
   return sqrt(deta*deta + dphi*dphi);
}

inline
double kinem::theta(double xyz1[3], double xyz2[3])
{
  double x=xyz2[0]-xyz1[0];
  double y=xyz2[1]-xyz1[1];
  double z=xyz2[2]-xyz1[2];
  return theta(x, y, z);
}

inline 
double kinem::theta(float xyz1[3], float xyz2[3]){
    double dxyz1[3]={xyz1[0],xyz1[1],xyz1[2]};
    double dxyz2[3]={xyz2[0],xyz2[1],xyz2[2]};
    return theta(dxyz1,dxyz2);
}

inline 
double kinem::theta(double x, double y, double z)
{
  return atan2(sqrt(x*x + y*y), z);
}

inline 
double kinem::theta(double etap)
{
  return 2.0*atan(exp(-etap));
}

inline
double kinem::eta(double xyz1[3], double xyz2[3])
{
  double x=xyz2[0]-xyz1[0];
  double y=xyz2[1]-xyz1[1];
  double z=xyz2[2]-xyz1[2];
  return eta(x, y, z);
}

inline
double kinem::eta(float xyz1[3], float xyz2[3]){
    double dxyz1[3]={xyz1[0],xyz1[1],xyz1[2]};
    double dxyz2[3]={xyz2[0],xyz2[1],xyz2[2]};
    return eta(dxyz1,dxyz2);
}

inline
double kinem::eta(double x, double y, double z)
{
  return 0.5*log((sqrt(x*x + y*y + z*z) + z + EPSILON) / 
                 (sqrt(x*x + y*y + z*z) - z + EPSILON));
}

inline
double kinem::eta(double th)
{
  if(th == 0) return ETA_LIMIT;
  if(th >= PI-0.0001) return -ETA_LIMIT;
  return -log(tan(th/2.0));
}

inline 
double kinem::y(double E, double pz)
{
  return 0.5 * log ((E+pz+EPSILON)/(E-pz+EPSILON));
}

inline
void kinem::uvectors(double u[3], double xyz1[3], double xyz2[3]){
  double xdiff=xyz2[0]-xyz1[0];
  double ydiff=xyz2[1]-xyz1[1];
  double zdiff=xyz2[2]-xyz1[2];
  double s=sqrt(xdiff*xdiff+ydiff*ydiff+zdiff*zdiff);
  if(s > 0) {
    u[0]=xdiff/s;
    u[1]=ydiff/s;
    u[2]=zdiff/s;
  }else{
    u[0]=0;
    u[1]=0;
    u[2]=0;
  }
}

inline 
void kinem::uvectors(float u[3], float xyz1[3], float xyz2[3]){
    double du[3];
    double dxyz1[3]={xyz1[0],xyz1[1],xyz1[2]};
    double dxyz2[3]={xyz2[0],xyz2[1],xyz2[2]};
    uvectors(du,dxyz1,dxyz2);
    u[0]=du[0];
    u[1]=du[1];
    u[2]=du[2];
}

inline
double kinem::tanl_from_theta(double theta){
  return tan(PI/2.0 - theta);
}

inline
double kinem::theta_from_tanl(double tanl){
  return PI/2.0 - atan(tanl);
}

#endif // INC_ANGLESUTIL
Revision:	1.1
Committed:	Tue Nov 11 23:01:21 2008 UTC (16 years, 5 months ago) by kukartse
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V00-03-01, ZMorph_BASE_20100408, gak040610_morphing, V00-02-02, gak011410, gak010310, ejterm2010_25nov2009, V00-02-01, V00-02-00, gak112409, CMSSW_22X_branch_base, segala101609, V00-01-15, V00-01-14, V00-01-13, V00-01-12, V00-01-11, V00-01-10, gak031009, gak030509, gak022309, gak021209, gak040209, gak012809, V00-01-09, V00-01-08, V00-01-07, V00-01-06, V00-01-05, V00-01-04, V00-00-07, V00-00-06, V00-00-05, V00-00-04, V00-01-03, V00-00-02, V00-00-01, HEAD
Branch point for:	ZMorph-V00-03-01, CMSSW_22X_branch
Log Message:	initial creation of a package for LJMET multivariate analysis
#	User	Rev	Content
1	kukartse	1.1	#ifndef INC_ANGLESUTIL
2			#define INC_ANGLESUTIL
3			///////////////////////////////////////////////////////////////////////////////
4			// File: AnglesUtil.hpp
5			//
6			// Purpose: Provide useful functions for calculating angles and eta
7			//
8			// Created: 4-NOV-1998 Serban Protopopescu
9			// History: replaced old KinemUtil
10			// Modified: 23-July-2000 Add rapidity calculation
11			// 14-Aug-2003 converted all functions to double precision
12			///////////////////////////////////////////////////////////////////////////////
13			// Dependencies (#includes)
14
15			#include <cmath>
16			#include <cstdlib>
17
18
19			namespace kinem{
20			const double PI=2.0*acos(0.);
21			const double TWOPI=2.0*PI;
22			const float ETA_LIMIT=15.0;
23			const float EPSILON=1.E-10;
24
25			// calculate phi from x, y
26			inline double phi(double x, double y);
27			// calculate phi for a line defined by xy1 and xy2 (xy2-xy1)
28			inline double phi(double xy1[2], double xy2[2]);
29			inline double phi(float xy1[2], float xy2[2]);
30
31			// calculate theta from x, y, z
32			inline double theta(double x, double y, double z);
33			// calculate theta for a line defined by xyz1 and xyz2 (xyz2-xyz1)
34			inline double theta(double xyz1[3], double xyz2[3]);
35			inline double theta(float xyz1[3], float xyz2[3]);
36			// calculate theta from eta
37			inline double theta(double etap);
38
39			// calculate eta from x, y, z (return also theta)
40			inline double eta(double x, double y, double z);
41			// calculate eta for a line defined by xyz1 and xyz2 (xyz2-xyz1)
42			inline double eta(double xyz1[3], double xyz2[3]);
43			inline double eta(float xyz1[3], float xyz2[3]);
44			// calculate eta from theta
45			inline double eta(double th);
46
47			// calculate rapidity from E, pz
48			inline double y(double E, double pz);
49
50			// calculate phi1-phi2 keeping value between 0 and pi
51			inline double delta_phi(double ph11, double phi2);
52			// calculate phi1-phi2 keeping value between -pi and pi
53			inline double signed_delta_phi(double ph11, double phi2);
54			// calculate deltaR
55			inline double delta_R(double eta1, double phi1, double eta2, double phi2);
56
57			// calculate unit vectors given two points
58			inline void uvectors(double u[3], double xyz1[3], double xyz2[3]);
59			inline void uvectors(float u[3], float xyz1[3], float xyz2[3]);
60
61			inline double tanl_from_theta(double theta);
62			inline double theta_from_tanl(double tanl);
63			}
64
65			inline
66			double kinem::phi(double x, double y)
67			{
68			double PHI=atan2(y, x);
69			return (PHI>=0)? PHI : kinem::TWOPI+PHI;
70			}
71
72			inline
73			double kinem::phi(float xy1[2], float xy2[2]){
74			double dxy1[2]={xy1[0],xy1[1]};
75			double dxy2[2]={xy2[0],xy2[1]};
76			return phi(dxy1,dxy2);
77			}
78
79			inline
80			double kinem::phi(double xy1[2], double xy2[2])
81			{
82			double x=xy2[0]-xy1[0];
83			double y=xy2[1]-xy1[1];
84			return phi(x, y);
85			}
86
87			inline
88			double kinem::delta_phi(double phi1, double phi2)
89			{
90			double PHI=fabs(phi1-phi2);
91			return (PHI<=PI)? PHI : kinem::TWOPI-PHI;
92			}
93
94			inline
95			double kinem::signed_delta_phi(double phi1, double phi2)
96			{
97			double phia=phi1;
98			if(phi1>PI) phia=phi1-kinem::TWOPI;
99			double phib=phi2;
100			if(phi2>PI) phib=phi2-kinem::TWOPI;
101			double dphi=phia-phib;
102			if(dphi>PI) dphi-=kinem::TWOPI;
103			if(dphi<-PI) dphi+=kinem::TWOPI;
104			return dphi;
105			}
106
107			inline double kinem::delta_R(double eta1, double phi1, double eta2, double phi2)
108			{
109			double deta = eta1-eta2;
110			double dphi = kinem::delta_phi(phi1,phi2);
111			return sqrt(detadeta + dphidphi);
112			}
113
114			inline
115			double kinem::theta(double xyz1[3], double xyz2[3])
116			{
117			double x=xyz2[0]-xyz1[0];
118			double y=xyz2[1]-xyz1[1];
119			double z=xyz2[2]-xyz1[2];
120			return theta(x, y, z);
121			}
122
123			inline
124			double kinem::theta(float xyz1[3], float xyz2[3]){
125			double dxyz1[3]={xyz1[0],xyz1[1],xyz1[2]};
126			double dxyz2[3]={xyz2[0],xyz2[1],xyz2[2]};
127			return theta(dxyz1,dxyz2);
128			}
129
130			inline
131			double kinem::theta(double x, double y, double z)
132			{
133			return atan2(sqrt(xx + yy), z);
134			}
135
136			inline
137			double kinem::theta(double etap)
138			{
139			return 2.0*atan(exp(-etap));
140			}
141
142			inline
143			double kinem::eta(double xyz1[3], double xyz2[3])
144			{
145			double x=xyz2[0]-xyz1[0];
146			double y=xyz2[1]-xyz1[1];
147			double z=xyz2[2]-xyz1[2];
148			return eta(x, y, z);
149			}
150
151			inline
152			double kinem::eta(float xyz1[3], float xyz2[3]){
153			double dxyz1[3]={xyz1[0],xyz1[1],xyz1[2]};
154			double dxyz2[3]={xyz2[0],xyz2[1],xyz2[2]};
155			return eta(dxyz1,dxyz2);
156			}
157
158			inline
159			double kinem::eta(double x, double y, double z)
160			{
161			return 0.5log((sqrt(xx + yy + zz) + z + EPSILON) /
162			(sqrt(xx + yy + z*z) - z + EPSILON));
163			}
164
165			inline
166			double kinem::eta(double th)
167			{
168			if(th == 0) return ETA_LIMIT;
169			if(th >= PI-0.0001) return -ETA_LIMIT;
170			return -log(tan(th/2.0));
171			}
172
173			inline
174			double kinem::y(double E, double pz)
175			{
176			return 0.5 * log ((E+pz+EPSILON)/(E-pz+EPSILON));
177			}
178
179			inline
180			void kinem::uvectors(double u[3], double xyz1[3], double xyz2[3]){
181			double xdiff=xyz2[0]-xyz1[0];
182			double ydiff=xyz2[1]-xyz1[1];
183			double zdiff=xyz2[2]-xyz1[2];
184			double s=sqrt(xdiffxdiff+ydiffydiff+zdiff*zdiff);
185			if(s > 0) {
186			u[0]=xdiff/s;
187			u[1]=ydiff/s;
188			u[2]=zdiff/s;
189			}else{
190			u[0]=0;
191			u[1]=0;
192			u[2]=0;
193			}
194			}
195
196			inline
197			void kinem::uvectors(float u[3], float xyz1[3], float xyz2[3]){
198			double du[3];
199			double dxyz1[3]={xyz1[0],xyz1[1],xyz1[2]};
200			double dxyz2[3]={xyz2[0],xyz2[1],xyz2[2]};
201			uvectors(du,dxyz1,dxyz2);
202			u[0]=du[0];
203			u[1]=du[1];
204			u[2]=du[2];
205			}
206
207			inline
208			double kinem::tanl_from_theta(double theta){
209			return tan(PI/2.0 - theta);
210			}
211
212			inline
213			double kinem::theta_from_tanl(double tanl){
214			return PI/2.0 - atan(tanl);
215			}
216
217			#endif // INC_ANGLESUTIL