ThreeVector.cc 8.36 KB
Newer Older
1
// -*- C++ -*-
2
// $Id: ThreeVector.cc,v 1.3 2003/08/13 20:00:14 garren Exp $
3
4
5
6
7
8
9
10
11
12
// ---------------------------------------------------------------------------
//
// This file is a part of the CLHEP - a Class Library for High Energy Physics.
//
// This is the implementation of the Hep3Vector class.
//
// See also ThreeVectorR.cc for implementation of Hep3Vector methods which 
// would couple in all the HepRotation methods.
//

13
#include "CLHEP/Vector/defs.h"
14
15
#include "CLHEP/Vector/ThreeVector.h"
#include "CLHEP/Vector/ZMxpv.h"
Lynn Garren's avatar
Lynn Garren committed
16
#include "CLHEP/Units/PhysicalConstants.h"
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43

#include <cmath>
#include <iostream>

namespace CLHEP  {

void Hep3Vector::setMag(double ma) {
  double factor = mag();
  if (factor == 0) {
    ZMthrowA ( ZMxpvZeroVector (
    "Hep3Vector::setMag : zero vector can't be stretched"));
  }else{
    factor = ma/factor;
    setX(x()*factor);
    setY(y()*factor);
    setZ(z()*factor);
  }
}

Hep3Vector & Hep3Vector::rotateUz(const Hep3Vector& NewUzVector) {
  // NewUzVector must be normalized !

  double u1 = NewUzVector.x();
  double u2 = NewUzVector.y();
  double u3 = NewUzVector.z();
  double up = u1*u1 + u2*u2;

44
45
46
47
48
49
50
51
52
53
54
  if (up > 0) {
    up = std::sqrt(up);
    double px = (u1 * u3 * x() - u2 * y()) / up + u1 * z();
    double py = (u2 * u3 * x() + u1 * y()) / up + u2 * z();
    double pz = -up * x() + u3 * z();
    set(px, py, pz);
  } else if (u3 < 0.) {
    setX(-x());
    setZ(-z());
  } // phi=0  teta=pi

55
56
57
58
  return *this;
}

double Hep3Vector::pseudoRapidity() const {
59
60
61
62
  double m1 = mag();
  if ( m1==  0   ) return  0.0;   
  if ( m1==  z() ) return  1.0E72;
  if ( m1== -z() ) return -1.0E72;
63
  return 0.5*std::log( (m1+z())/(m1-z()) );
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
}

std::ostream & operator<< (std::ostream & os, const Hep3Vector & v) {
  return os << "(" << v.x() << "," << v.y() << "," << v.z() << ")";
}

void ZMinput3doubles ( std::istream & is, const char * type,
                       double & x, double & y, double & z );

std::istream & operator>>(std::istream & is, Hep3Vector & v) {
  double x, y, z;
  ZMinput3doubles ( is, "Hep3Vector", x, y, z );
  v.set(x, y, z);
  return  is;
}  // operator>>()

const Hep3Vector HepXHat(1.0, 0.0, 0.0);
const Hep3Vector HepYHat(0.0, 1.0, 0.0);
const Hep3Vector HepZHat(0.0, 0.0, 1.0);

//-------------------
//
// New methods introduced when ZOOM PhysicsVectors was merged in:
//
//-------------------

90
Hep3Vector & Hep3Vector::rotateX (double phi1) {
91
92
  double sinphi = std::sin(phi1);
  double cosphi = std::cos(phi1);
93
94
95
96
  double ty = y() * cosphi - z() * sinphi;
  double tz = z() * cosphi + y() * sinphi;
  setY(ty);
  setZ(tz);
97
98
99
  return *this;
} /* rotateX */

100
Hep3Vector & Hep3Vector::rotateY (double phi1) {
101
102
  double sinphi = std::sin(phi1);
  double cosphi = std::cos(phi1);
103
104
105
106
  double tx = x() * cosphi + z() * sinphi;
  double tz = z() * cosphi - x() * sinphi;
  setX(tx);
  setZ(tz);
107
108
109
  return *this;
} /* rotateY */

110
Hep3Vector & Hep3Vector::rotateZ (double phi1) {
111
112
  double sinphi = std::sin(phi1);
  double cosphi = std::cos(phi1);
113
114
115
116
  double tx = x() * cosphi - y() * sinphi;
  double ty = y() * cosphi + x() * sinphi;
  setX(tx);
  setY(ty);
117
118
119
120
121
122
123
124
125
126
127
128
129
  return *this;
} /* rotateZ */

bool Hep3Vector::isNear(const Hep3Vector & v, double epsilon) const {
  double limit = dot(v)*epsilon*epsilon;
  return ( (*this - v).mag2() <= limit );
} /* isNear() */

double Hep3Vector::howNear(const Hep3Vector & v ) const {
  // | V1 - V2 | **2  / V1 dot V2, up to 1
  double d   = (*this - v).mag2();
  double vdv = dot(v);
  if ( (vdv > 0) && (d < vdv)  ) {
130
    return std::sqrt (d/vdv);
131
132
133
134
135
136
137
138
139
  } else if ( (vdv == 0) && (d == 0) ) {
    return 0;
  } else {
    return 1;
  }
} /* howNear */

double Hep3Vector::deltaPhi  (const Hep3Vector & v2) const {
  double dphi = v2.getPhi() - getPhi();
Lynn Garren's avatar
Lynn Garren committed
140
141
142
143
  if ( dphi > CLHEP::pi ) {
    dphi -= CLHEP::twopi;
  } else if ( dphi <= -CLHEP::pi ) {
    dphi += CLHEP::twopi;
144
145
146
147
148
149
150
  }
  return dphi;
} /* deltaPhi */

double Hep3Vector::deltaR ( const Hep3Vector & v ) const {
  double a = eta() - v.eta();
  double b = deltaPhi(v); 
151
  return std::sqrt ( a*a + b*b );
152
153
154
155
156
157
158
159
} /* deltaR */

double Hep3Vector::cosTheta(const Hep3Vector & q) const {
  double arg;
  double ptot2 = mag2()*q.mag2();
  if(ptot2 <= 0) {
    arg = 0.0;
  }else{
160
    arg = dot(q)/std::sqrt(ptot2);
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
    if(arg >  1.0) arg =  1.0;
    if(arg < -1.0) arg = -1.0;
  }
  return arg;
}

double Hep3Vector::cos2Theta(const Hep3Vector & q) const {
  double arg;
  double ptot2 = mag2();
  double qtot2 = q.mag2();
  if ( ptot2 == 0 || qtot2 == 0 )  {
    arg = 1.0;
  }else{
    double pdq = dot(q);
    arg = (pdq/ptot2) * (pdq/qtot2);
        // More naive methods overflow on vectors which can be squared
        // but can't be raised to the 4th power.
    if(arg >  1.0) arg =  1.0;
 }
 return arg;
}

183
184
185
void Hep3Vector::setEta (double eta1) {
  double phi1 = 0;
  double r1;
186
187
  if ( (x() == 0) && (y() == 0) ) {
    if (z() == 0) {
188
189
190
191
192
193
      ZMthrowC (ZMxpvZeroVector(
        "Attempt to set eta of zero vector -- vector is unchanged"));
      return;
    }
    ZMthrowC (ZMxpvZeroVector(
      "Attempt to set eta of vector along Z axis -- will use phi = 0"));
194
    r1 = std::fabs(z());
195
  } else {
196
197
    r1 = getR();
    phi1 = getPhi();
198
  }
199
  double tanHalfTheta = std::exp ( -eta1 );
200
  double cosTheta1 =
201
        (1 - tanHalfTheta*tanHalfTheta) / (1 + tanHalfTheta*tanHalfTheta);
202
  double rho1 = r1*std::sqrt(1 - cosTheta1*cosTheta1);
203
204
205
  setZ(r1 * cosTheta1);
  setY(rho1 * std::sin (phi1));
  setX(rho1 * std::cos (phi1));
206
207
208
  return;
}

209
void Hep3Vector::setCylTheta (double theta1) {
210
211
212

  // In cylindrical coords, set theta while keeping rho and phi fixed

213
214
  if ( (x() == 0) && (y() == 0) ) {
    if (z() == 0) {
215
216
217
218
      ZMthrowC (ZMxpvZeroVector(
        "Attempt to set cylTheta of zero vector -- vector is unchanged"));
      return;
    }
219
    if (theta1 == 0) {
220
      setZ(std::fabs(z()));
221
222
      return;
    }
223
    if (theta1 == CLHEP::pi) {
224
      setZ(-std::fabs(z()));
225
226
227
228
229
230
      return;
    }
    ZMthrowC (ZMxpvZeroVector(
      "Attempt set cylindrical theta of vector along Z axis "
      "to a non-trivial value, while keeping rho fixed -- "
      "will return zero vector"));
231
    setZ(0.0);
232
233
    return;
  }
234
  if ( (theta1 < 0) || (theta1 > CLHEP::pi) ) {
235
236
237
238
    ZMthrowC (ZMxpvUnusualTheta(
      "Setting Cyl theta of a vector based on a value not in [0, PI]"));
        // No special return needed if warning is ignored.
  }
239
240
241
  double phi1 (getPhi());
  double rho1 = getRho();
  if ( (theta1 == 0) || (theta1 == CLHEP::pi) ) {
242
243
244
    ZMthrowC (ZMxpvInfiniteVector(
      "Attempt to set cylindrical theta to 0 or PI "
      "while keeping rho fixed -- infinite Z will be computed"));
245
      setZ((theta1==0) ? 1.0E72 : -1.0E72);
246
247
    return;
  }
248
249
250
  setZ(rho1 / std::tan (theta1));
  setY(rho1 * std::sin (phi1));
  setX(rho1 * std::cos (phi1));
251
252
253

} /* setCylTheta */

254
void Hep3Vector::setCylEta (double eta1) {
255
256
257

  // In cylindrical coords, set eta while keeping rho and phi fixed

258
  double theta1 = 2 * std::atan ( std::exp (-eta1) );
259
260
261
262
263
264

        //-| The remaining code is similar to setCylTheta,  The reason for
        //-| using a copy is so as to be able to change the messages in the
        //-| ZMthrows to say eta rather than theta.  Besides, we assumedly
        //-| need not check for theta of 0 or PI.

265
266
  if ( (x() == 0) && (y() == 0) ) {
    if (z() == 0) {
267
268
269
270
      ZMthrowC (ZMxpvZeroVector(
        "Attempt to set cylEta of zero vector -- vector is unchanged"));
      return;
    }
271
    if (theta1 == 0) {
272
      setZ(std::fabs(z()));
273
274
      return;
    }
275
    if (theta1 == CLHEP::pi) {
276
      setZ(-std::fabs(z()));
277
278
279
280
281
282
      return;
    }
    ZMthrowC (ZMxpvZeroVector(
      "Attempt set cylindrical eta of vector along Z axis "
      "to a non-trivial value, while keeping rho fixed -- "
      "will return zero vector"));
283
    setZ(0.0);
284
285
    return;
  }
286
287
  double phi1 (getPhi());
  double rho1 = getRho();
288
289
290
  setZ(rho1 / std::tan (theta1));
  setY(rho1 * std::sin (phi1));
  setX(rho1 * std::cos (phi1));
291
292
293
294

} /* setCylEta */


295
Hep3Vector operator/ ( const Hep3Vector & v1, double c ) {
296
297
298
299
300
  if (c == 0) {
    ZMthrowA ( ZMxpvInfiniteVector (
      "Attempt to divide vector by 0 -- "
      "will produce infinities and/or NANs"));
  } 
301
  return v1 * (1.0/c);
302
303
304
305
306
307
308
309
} /* v / c */

Hep3Vector & Hep3Vector::operator/= (double c) {
  if (c == 0) {
    ZMthrowA (ZMxpvInfiniteVector(
      "Attempt to do vector /= 0 -- "
      "division by zero would produce infinite or NAN components"));
  }
310
  *this *= 1.0/c;
311
312
313
314
315
316
317
318
319
  return *this;
}

double Hep3Vector::tolerance = Hep3Vector::ToleranceTicks * 2.22045e-16;

}  // namespace CLHEP