Logo ROOT   6.10/09
Reference Guide
coordinates3D.cxx
Go to the documentation of this file.
1 // $Id $
2 //
3 // Tests that each form of vector has all the properties that stem from
4 // owning and forwarding to a coordinates instance, and that they give proper
5 // results.
6 //
7 // 3D vectors have:
8 //
9 // accessors x(), y(), z(), r(), theta(), phi(), rho(), eta()
10 //
11 // =================================================================
12 
16 #include "Math/GenVector/Polar3D.h"
19 #include "Math/GenVector/etaMax.h"
20 
21 #include "Math/Vector3Dfwd.h"
22 
23 #include "CoordinateTraits.h"
24 
25 #include <iostream>
26 #include <limits>
27 #include <cmath>
28 #include <vector>
29 
30 //#define TRACE1
31 #define DEBUG
32 
33 using namespace ROOT::Math;
34 
35 
36 
37 
38 template <typename T1, typename T2 >
39 struct Precision {
40  enum { result = std::numeric_limits<T1>::digits <= std::numeric_limits<T2>::digits };
41 };
42 
43 template <typename T1, typename T2, bool>
44 struct LessPreciseType {
45  typedef T1 type;
46 };
47 template <typename T1, typename T2>
48 struct LessPreciseType<T1, T2, false> {
49  typedef T2 type;
50 };
51 
52 
53 template <typename Scalar1, typename Scalar2>
54 int
55 closeEnough ( Scalar1 s1, Scalar2 s2, std::string const & coord, double ticks ) {
56  int ret = 0;
57  int pr = std::cout.precision(18);
61  Scalar epsilon = (eps1 >= eps2) ? eps1 : eps2;
62  Scalar ss1 (s1);
63  Scalar ss2 (s2);
64  Scalar diff = ss1 - ss2;
65  if (diff < 0) diff = -diff;
66  if (ss1 == 0 || ss2 == 0) { // TODO - the ss2==0 makes a big change??
67  if ( diff > ticks*epsilon ) {
68  ret=3;
69  std::cout << "\nAbsolute discrepancy in " << coord << "(): "
70  << ss1 << " != " << ss2 << "\n"
71  << " (Allowed discrepancy is " << ticks*epsilon
72  << ")\nDifference is " << diff/epsilon << " ticks\n";
73  }
74  std::cout.precision (pr);
75  return ret;
76  }
77  // infinity dicrepancy musy be checked with max precision
78  long double sd1(ss1);
79  long double sd2(ss2);
80  if ( (sd1 + sd2 == sd1) != (sd1 + sd2 == sd2) ) {
81  ret=5;
82  std::cout << "\nInfinity discrepancy in " << coord << "(): "
83  << sd1 << " != " << sd2 << "\n";
84  std::cout.precision (pr);
85  return ret;
86  }
87  Scalar denom = ss1 > 0 ? ss1 : -ss1;
88  if ((diff/denom > ticks*epsilon) && (diff > ticks*epsilon)) {
89  ret=9;
90  std::cout << "\nDiscrepancy in " << coord << "(): "
91  << ss1 << " != " << ss2 << "\n"
92  << " (Allowed discrepancy is " << ticks*epsilon*ss1
93  << ")\nDifference is " << (diff/denom)/epsilon << " ticks\n";
94  }
95  std::cout.precision (pr);
96  return ret;
97 }
98 
99 template <class V1, class V2>
100 int compare3D (const V1 & v1, const V2 & v2, double ticks) {
101  int ret =0;
102  typedef typename V1::CoordinateType CoordType1;
103  typedef typename V2::CoordinateType CoordType2;
104  ret |= closeEnough ( v1.x(), v2.x(), "x" ,ticks);
105  ret |= closeEnough ( v1.y(), v2.y(), "y" ,ticks);
106  ret |= closeEnough ( v1.z(), v2.z(), "z" ,ticks);
107  ret |= closeEnough ( v1.rho(), v2.rho(), "rho" ,ticks);
108  // case in phi that difference is close to 2pi
109  typedef typename V2::Scalar Scalar;
110  Scalar phi2 = v2.phi();
111  if (std::abs(v1.phi()- phi2 ) > ROOT::Math::Pi() ) {
112  if (phi2<0)
113  phi2 += 2.*ROOT::Math::Pi();
114  else
115  phi2 -= 2*ROOT::Math::Pi();
116  }
117  ret |= closeEnough ( v1.phi(), phi2, "phi" ,ticks);
118  ret |= closeEnough ( v1.r(), v2.r(), "r" ,ticks);
119  ret |= closeEnough ( v1.theta(), v2.theta(), "theta" ,ticks);
120  ret |= closeEnough ( v1.mag2(), v2.mag2(), "mag2" ,ticks);
121  ret |= closeEnough ( v1.perp2(), v2.perp2(), "perp2" ,ticks);
122  if ( v1.rho() > 0 && v2.rho() > 0 ) { // eta can legitimately vary if rho == 0
123  ret |= closeEnough ( v1.eta(), v2.eta(), "eta" ,ticks);
124  }
125 
126  if (ret != 0) {
127  std::cout << "\nDiscrepancy detected (see above) is between:\n "
128  << CoordinateTraits<CoordType1>::name() << " and\n "
130  << "with v = (" << v1.x() << ", " << v1.y() << ", " << v1.z()
131  << ")\n\n\n";
132  }
133  else {
134 #ifdef DEBUG
135  std::cout << ".";
136 #endif
137  }
138 
139 
140  return ret;
141 }
142 
143 template <class C>
144 int test3D ( const DisplacementVector3D<C> & v, double ticks ) {
145 
146 #ifdef DEBUG
147  std::cout <<"\n>>>>> Testing 3D from " << v << " ticks = " << ticks << std::endl;
148 #endif
149 
150  int ret = 0;
151  DisplacementVector3D< Cartesian3D<double> > vxyz_d (v.x(), v.y(), v.z());
152 
153  double r = std::sqrt (v.x()*v.x() + v.y()*v.y() + v.z()*v.z());
154  double rho = std::sqrt (v.x()*v.x() + v.y()*v.y());
155  double z = v.z();
156 // double theta = r>0 ? std::acos ( z/r ) : 0;
157 // if (std::abs( std::abs(z) - r) < 10*r* std::numeric_limits<double>::epsilon() )
158  double theta = std::atan2( rho, z ); // better when theta is small or close to pi
159 
160  double phi = rho>0 ? std::atan2 (v.y(), v.x()) : 0;
161  DisplacementVector3D< Polar3D<double> > vrtp_d ( r, theta, phi );
162 
163  double eta;
164  if (rho != 0) {
165  eta = -std::log(std::tan(theta/2));
166  #ifdef TRACE1
167  std::cout << ":::: rho != 0\n"
168  << ":::: theta = << " << theta
169  <<"/n:::: tan(theta/2) = " << std::tan(theta/2)
170  <<"\n:::: eta = " << eta << "\n";
171  #endif
172  } else if (v.z() == 0) {
173  eta = 0;
174  #ifdef TRACE1
175  std::cout << ":::: v.z() == 0\n"
176  <<"\n:::: eta = " << eta << "\n";
177  #endif
178  } else if (v.z() > 0) {
179  eta = v.z() + etaMax<long double>();
180  #ifdef TRACE1
181  std::cout << ":::: v.z() > 0\n"
182  << ":::: etaMax = " << etaMax<long double>()
183  <<"\n:::: eta = " << eta << "\n";
184  #endif
185  } else {
186  eta = v.z() - etaMax<long double>();
187  #ifdef TRACE1
188  std::cout << ":::: v.z() < 0\n"
189  << ":::: etaMax = " << etaMax<long double>()
190  <<"\n:::: eta = " << eta << "\n";
191  #endif
192  }
193 
194 #ifdef DEBUG
195  std::cout << " Testing DisplacementVector3D : ";
196 #endif
197 
198  DisplacementVector3D< CylindricalEta3D<double> > vrep_d ( rho, eta, phi );
199 
200  ret |= compare3D( vxyz_d, vrtp_d, ticks);
201  ret |= compare3D( vrtp_d, vxyz_d, ticks);
202  ret |= compare3D( vxyz_d, vrep_d, ticks);
203  ret |= compare3D( vrtp_d, vrep_d, ticks);
204 
205  DisplacementVector3D< Cylindrical3D<double> > vrzp_d ( rho, z, phi );
206 
207  ret |= compare3D( vxyz_d, vrzp_d, ticks);
208  ret |= compare3D( vrtp_d, vrzp_d, ticks);
209  ret |= compare3D( vrep_d, vrzp_d, ticks);
210 
211  DisplacementVector3D< Cartesian3D<float> > vxyz_f (v.x(), v.y(), v.z());
212  DisplacementVector3D< Polar3D<float> > vrtp_f ( r, theta, phi );
213  DisplacementVector3D< CylindricalEta3D<float> > vrep_f ( rho, eta, phi );
214 
215  ret |= compare3D( vxyz_d, vxyz_f, ticks);
216  ret |= compare3D( vxyz_d, vrep_f, ticks);
217  ret |= compare3D( vrtp_d, vrep_f, ticks);
218 
219  ret |= compare3D( vxyz_f, vxyz_f, ticks);
220  ret |= compare3D( vxyz_f, vrep_f, ticks);
221  ret |= compare3D( vrtp_f, vrep_f, ticks);
222 
223 #ifdef DEBUG
224  if (ret == 0) std::cout << "\t OK\n";
225  else {
226  std::cout << "\t FAIL\n";
227  std::cerr << "\n>>>>> Testing DisplacementVector3D from " << v << " ticks = " << ticks
228  << "\t:\t FAILED\n";
229  }
230  std::cout << " Testing PositionVector3D : ";
231 #endif
232 
233 
234  PositionVector3D< Cartesian3D <double> > pxyz_d; pxyz_d = vxyz_d;
235  PositionVector3D< Polar3D <double> > prtp_d; prtp_d = vrtp_d;
236  PositionVector3D< CylindricalEta3D<double> > prep_d; prep_d = vrep_d;
237  PositionVector3D< Cylindrical3D <double> > przp_d; przp_d = vrzp_d;
238 
239  ret |= compare3D( pxyz_d, prtp_d, ticks);
240  ret |= compare3D( vxyz_d, prep_d, ticks);
241  ret |= compare3D( vrtp_d, prep_d, ticks);
242  ret |= compare3D( vxyz_d, przp_d, ticks);
243 
244  PositionVector3D< Cartesian3D<float> > pxyz_f (v.x(), v.y(), v.z());
245  PositionVector3D< Polar3D<float> > prtp_f ( r, theta, phi );
246  PositionVector3D< CylindricalEta3D<float> > prep_f ( rho, eta, phi );
247 
248  ret |= compare3D( vxyz_d, pxyz_f, ticks);
249  ret |= compare3D( vxyz_d, prep_f, ticks);
250  ret |= compare3D( vrtp_d, prep_f, ticks);
251 
252  ret |= compare3D( vxyz_f, pxyz_f, ticks);
253  ret |= compare3D( vxyz_f, prep_f, ticks);
254  ret |= compare3D( vrtp_f, prep_f, ticks);
255 
256 #ifdef DEBUG
257  if (ret == 0) std::cout << "\t\t OK\n";
258  else {
259  std::cout << "\t FAIL\n";
260  std::cerr << "\n>>>>> Testing PositionVector3D from " << v << " ticks = " << ticks
261  << "\t:\t FAILED\n";
262  }
263 #endif
264  return ret;
265 }
266 
267 
268 
270  int ret = 0;
271 
272  ret |= test3D (XYZVector ( 0.0, 0.0, 0.0 ) ,2 );
273  ret |= test3D (XYZVector ( 1.0, 2.0, 3.0 ) ,6 );
274  ret |= test3D (XYZVector ( -1.0, 2.0, 3.0 ) ,6 );
275  ret |= test3D (XYZVector ( 1.0, -2.0, 3.0 ) ,6 );
276  ret |= test3D (XYZVector ( 1.0, 2.0, -3.0 ) ,6 );
277  ret |= test3D (XYZVector ( -1.0, -2.0, 3.0 ) ,6 );
278  ret |= test3D (XYZVector ( -1.0, 2.0, -3.0 ) ,6 );
279  ret |= test3D (XYZVector ( 1.0, -2.0, -3.0 ) ,6 );
280  ret |= test3D (XYZVector ( -1.0, -2.0, -3.0 ) ,6 );
281  ret |= test3D (XYZVector ( 8.0, 0.0, 0.0 ) ,6 );
282  ret |= test3D (XYZVector ( -8.0, 0.0, 0.0 ) ,12 );
283  ret |= test3D (XYZVector ( 0.0, 9.0, 0.0 ) ,6 );
284  ret |= test3D (XYZVector ( 0.0, -9.0, 0.0 ) ,6 );
285 // rho == 0 tests the beyon-eta-max cases of cylindricalEta
286  ret |= test3D (XYZVector ( 0.0, 0.0, 7.0 ) ,8 );
287  ret |= test3D (XYZVector ( 0.0, 0.0, -7.0 ) ,8 );
288 // Larger ratios among coordinates presents a precision challenge
289  ret |= test3D (XYZVector ( 16.0, 0.02, .01 ) ,10 );
290  ret |= test3D (XYZVector ( -16.0, 0.02, .01 ) ,10 );
291  ret |= test3D (XYZVector ( -.01, 16.0, .01 ) ,2000 );
292  ret |= test3D (XYZVector ( -.01, -16.0, .01 ) ,2000 );
293  ret |= test3D (XYZVector ( 1.0, 2.0, 30.0 ) ,10 );
294  // NOTE -- these larger erros are likely the results of treating
295  // the vector in a ctor or assignment as foreign...
296  // NO -- I'm fouling up the value of x() !!!!!
297 // As we push to higher z with zero rho, some accuracy loss is expected
298  ret |= test3D (XYZVector ( 0.0, 0.0, 15.0 ) ,30 );
299  ret |= test3D (XYZVector ( 0.0, 0.0, -15.0 ) ,30 );
300  ret |= test3D (XYZVector ( 0.0, 0.0, 35.0 ) ,30 );
301  ret |= test3D (XYZVector ( 0.0, 0.0, -35.0 ) ,30 );
302 // When z is big compared to rho, it is very hard to get precision in polar/eta:
303  ret |= test3D (XYZVector ( 0.01, 0.02, 16.0 ) ,10 );
304  ret |= test3D (XYZVector ( 0.01, 0.02, -16.0 ) ,40000 );
305 // test case when eta is large
306  ret |= test3D (XYZVector ( 1.E-8, 1.E-8, 10.0 ) , 20 );
307 // when z is neg error is larger in eta when calculated from polar
308 // since we have a larger error in theta which is closer to pi
309  ret |= test3D (XYZVector ( 1.E-8, 1.E-8, -10.0 ) ,2.E9 );
310 
311  // small value of z
312  ret |= test3D (XYZVector ( 10., 10., 1.E-8 ) ,1.0E6 );
313  ret |= test3D (XYZVector ( 10., 10., -1.E-8 ) ,1.0E6 );
314 
315 
316  return ret;
317 }
318 
319 int main() {
320  int ret = coordinates3D();
321  if (ret) std::cerr << "test FAILED !!! " << std::endl;
322  else std::cout << "test OK " << std::endl;
323  return ret;
324 }
int coordinates3D()
Class describing a generic position vector (point) in 3 dimensions.
int test3D(const DisplacementVector3D< C > &v, double ticks)
double sqrt(double)
static const std::string name()
int compare3D(const V1 &v1, const V2 &v2, double ticks)
Class describing a generic displacement vector in 3 dimensions.
TRandom2 r(17)
SVector< double, 2 > v
Definition: Dict.h:5
double Pi()
Mathematical constants.
Definition: Math.h:90
int closeEnough(Scalar1 s1, Scalar2 s2, std::string const &coord, double ticks)
REAL epsilon
Definition: triangle.c:617
constexpr Double_t E()
Definition: TMath.h:74
RooCmdArg Precision(Double_t prec)
double atan2(double, double)
int type
Definition: TGX11.cxx:120
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
#define T1
Definition: md5.inl:145
double tan(double)
int main()
double result[121]
Rotation3D::Scalar Scalar
double log(double)