Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/externals/clhep/src/LorentzVector.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /externals/clhep/src/LorentzVector.cc (Version 11.3.0) and /externals/clhep/src/LorentzVector.cc (Version 11.1.3)


  1 // -*- C++ -*-                                      1 // -*- C++ -*-
  2 // -------------------------------------------      2 // ---------------------------------------------------------------------------
  3 //                                                  3 //
  4 // This file is a part of the CLHEP - a Class       4 // This file is a part of the CLHEP - a Class Library for High Energy Physics.
  5 //                                                  5 //
  6 // This is the implementation of that portion       6 // This is the implementation of that portion of the HepLorentzVector class
  7 // which was in the original CLHEP and which d      7 // which was in the original CLHEP and which does not force loading of either
  8 // Rotation.cc or LorentzRotation.cc                8 // Rotation.cc or LorentzRotation.cc
  9 //                                                  9 //
 10                                                    10 
 11 #include "CLHEP/Vector/LorentzVector.h"            11 #include "CLHEP/Vector/LorentzVector.h"
 12                                                    12 
 13 #include <cmath>                                   13 #include <cmath>
 14 #include <iostream>                                14 #include <iostream>
 15                                                    15 
 16 namespace CLHEP  {                                 16 namespace CLHEP  {
 17                                                    17 
 18 double HepLorentzVector::tolerance =               18 double HepLorentzVector::tolerance = 
 19         Hep3Vector::ToleranceTicks * 2.22045e-     19         Hep3Vector::ToleranceTicks * 2.22045e-16;
 20 double HepLorentzVector::metric = 1.0;             20 double HepLorentzVector::metric = 1.0;
 21                                                    21 
 22 double HepLorentzVector::operator () (int i) c     22 double HepLorentzVector::operator () (int i) const {
 23   switch(i) {                                      23   switch(i) {
 24   case X:                                          24   case X:
 25   case Y:                                          25   case Y:
 26   case Z:                                          26   case Z:
 27     return pp(i);                                  27     return pp(i);
 28   case T:                                          28   case T:
 29     return e();                                    29     return e();
 30   default:                                         30   default:
 31     std::cerr << "HepLorentzVector subscriptin     31     std::cerr << "HepLorentzVector subscripting: bad index (" << i << ")"
 32      << std::endl;                                 32      << std::endl;
 33   }                                                33   }
 34   return 0.;                                       34   return 0.;
 35 }                                                  35 }  
 36                                                    36 
 37 double & HepLorentzVector::operator () (int i)     37 double & HepLorentzVector::operator () (int i) {
 38   static double dummy;                             38   static double dummy;
 39   switch(i) {                                      39   switch(i) {
 40   case X:                                          40   case X:
 41   case Y:                                          41   case Y:
 42   case Z:                                          42   case Z:
 43     return pp(i);                                  43     return pp(i);
 44   case T:                                          44   case T:
 45     return ee;                                     45     return ee;
 46   default:                                         46   default:
 47     std::cerr                                      47     std::cerr
 48       << "HepLorentzVector subscripting: bad i     48       << "HepLorentzVector subscripting: bad index (" << i << ")"
 49       << std::endl;                                49       << std::endl;
 50     return dummy;                                  50     return dummy;
 51   }                                                51   }
 52 }                                                  52 }
 53                                                    53 
 54 HepLorentzVector & HepLorentzVector::boost         54 HepLorentzVector & HepLorentzVector::boost
 55         (double bx, double by, double bz){         55         (double bx, double by, double bz){
 56   double b2 = bx*bx + by*by + bz*bz;               56   double b2 = bx*bx + by*by + bz*bz;
 57   double ggamma = 1.0 / std::sqrt(1.0 - b2);       57   double ggamma = 1.0 / std::sqrt(1.0 - b2);
 58   double bp = bx*x() + by*y() + bz*z();            58   double bp = bx*x() + by*y() + bz*z();
 59   double gamma2 = b2 > 0 ? (ggamma - 1.0)/b2 :     59   double gamma2 = b2 > 0 ? (ggamma - 1.0)/b2 : 0.0;
 60                                                    60 
 61   setX(x() + gamma2*bp*bx + ggamma*bx*t());        61   setX(x() + gamma2*bp*bx + ggamma*bx*t());
 62   setY(y() + gamma2*bp*by + ggamma*by*t());        62   setY(y() + gamma2*bp*by + ggamma*by*t());
 63   setZ(z() + gamma2*bp*bz + ggamma*bz*t());        63   setZ(z() + gamma2*bp*bz + ggamma*bz*t());
 64   setT(ggamma*(t() + bp));                         64   setT(ggamma*(t() + bp));
 65   return *this;                                    65   return *this;
 66 }                                                  66 }
 67                                                    67 
 68 HepLorentzVector & HepLorentzVector::rotateX(d     68 HepLorentzVector & HepLorentzVector::rotateX(double a) {
 69   pp.rotateX(a);                                   69   pp.rotateX(a); 
 70   return *this;                                    70   return *this; 
 71 }                                                  71 }
 72 HepLorentzVector & HepLorentzVector::rotateY(d     72 HepLorentzVector & HepLorentzVector::rotateY(double a) { 
 73   pp.rotateY(a);                                   73   pp.rotateY(a); 
 74   return *this;                                    74   return *this; 
 75 }                                                  75 }
 76 HepLorentzVector & HepLorentzVector::rotateZ(d     76 HepLorentzVector & HepLorentzVector::rotateZ(double a) { 
 77   pp.rotateZ(a);                                   77   pp.rotateZ(a); 
 78   return *this;                                    78   return *this; 
 79 }                                                  79 }
 80                                                    80 
 81 HepLorentzVector & HepLorentzVector::rotateUz(     81 HepLorentzVector & HepLorentzVector::rotateUz(const Hep3Vector &v1) {
 82   pp.rotateUz(v1);                                 82   pp.rotateUz(v1);
 83   return *this;                                    83   return *this;
 84 }                                                  84 }
 85                                                    85 
 86 std::ostream & operator<< (std::ostream & os,      86 std::ostream & operator<< (std::ostream & os, const HepLorentzVector & v1)
 87 {                                                  87 {
 88   return os << "(" << v1.x() << "," << v1.y()      88   return os << "(" << v1.x() << "," << v1.y() << "," << v1.z()
 89       << ";" << v1.t() << ")";                     89       << ";" << v1.t() << ")";
 90 }                                                  90 }
 91                                                    91 
 92 std::istream & operator>> (std::istream & is,      92 std::istream & operator>> (std::istream & is, HepLorentzVector & v1) {
 93                                                    93 
 94 // Required format is ( a, b, c; d ) that is,      94 // Required format is ( a, b, c; d ) that is, four numbers, preceded by
 95 // (, followed by ), components of the spatial     95 // (, followed by ), components of the spatial vector separated by commas,
 96 // time component separated by semicolon. The      96 // time component separated by semicolon. The four numbers are taken
 97 // as x, y, z, t.                                  97 // as x, y, z, t.
 98                                                    98 
 99   double x, y, z, t;                               99   double x, y, z, t;
100   char c;                                         100   char c;
101                                                   101 
102   is >> std::ws >> c;                             102   is >> std::ws >> c;
103     // ws is defined to invoke eatwhite(istrea    103     // ws is defined to invoke eatwhite(istream & )
104     // see (Stroustrup gray book) page 333 and    104     // see (Stroustrup gray book) page 333 and 345.
105   if (is.fail() || c != '(' ) {                   105   if (is.fail() || c != '(' ) {
106     std::cerr << "Could not find required open    106     std::cerr << "Could not find required opening parenthesis "
107         << "in input of a HepLorentzVector" <<    107         << "in input of a HepLorentzVector" << std::endl;
108     return is;                                    108     return is;
109   }                                               109   }
110                                                   110 
111   is >> x >> std::ws >> c;                        111   is >> x >> std::ws >> c;
112   if (is.fail() || c != ',' ) {                   112   if (is.fail() || c != ',' ) {
113     std::cerr << "Could not find x value and r    113     std::cerr << "Could not find x value and required trailing comma "
114         << "in input of a HepLorentzVector" <<    114         << "in input of a HepLorentzVector" << std::endl; 
115     return is;                                    115     return is;
116   }                                               116   }
117                                                   117 
118   is >> y >> std::ws >> c;                        118   is >> y >> std::ws >> c;
119   if (is.fail() || c != ',' ) {                   119   if (is.fail() || c != ',' ) {
120     std::cerr << "Could not find y value and r    120     std::cerr << "Could not find y value and required trailing comma "
121               <<  "in input of a HepLorentzVec    121               <<  "in input of a HepLorentzVector" << std::endl;
122     return is;                                    122     return is;
123   }                                               123   }
124                                                   124 
125   is >> z >> std::ws >> c;                        125   is >> z >> std::ws >> c;
126   if (is.fail() || c != ';' ) {                   126   if (is.fail() || c != ';' ) {
127     std::cerr << "Could not find z value and r    127     std::cerr << "Could not find z value and required trailing semicolon "
128      <<  "in input of a HepLorentzVector" << s    128      <<  "in input of a HepLorentzVector" << std::endl;
129     return is;                                    129     return is;
130   }                                               130   }
131                                                   131 
132   is >> t >> std::ws >> c;                        132   is >> t >> std::ws >> c;
133   if (is.fail() || c != ')' ) {                   133   if (is.fail() || c != ')' ) {
134     std::cerr << "Could not find t value and r    134     std::cerr << "Could not find t value and required close parenthesis "
135      << "in input of a HepLorentzVector" << st    135      << "in input of a HepLorentzVector" << std::endl;
136     return is;                                    136     return is;
137   }                                               137   }
138                                                   138 
139   v1.setX(x);                                     139   v1.setX(x);
140   v1.setY(y);                                     140   v1.setY(y);
141   v1.setZ(z);                                     141   v1.setZ(z);
142   v1.setT(t);                                     142   v1.setT(t);
143   return is;                                      143   return is;
144 }                                                 144 }
145                                                   145 
146 // The following were added when ZOOM classes     146 // The following were added when ZOOM classes were merged in:
147                                                   147 
148 HepLorentzVector & HepLorentzVector::operator     148 HepLorentzVector & HepLorentzVector::operator /= (double c) {
149 //  if (c == 0) {                                 149 //  if (c == 0) {
150 //    std::cerr << "HepLorentzVector::operator    150 //    std::cerr << "HepLorentzVector::operator /=() - "
151 //      << "Attempt to do LorentzVector /= 0 -    151 //      << "Attempt to do LorentzVector /= 0 -- \n"
152 //      << "division by zero would produce inf    152 //      << "division by zero would produce infinite or NAN components"
153 //      << std::endl;                             153 //      << std::endl;
154 //  }                                             154 //  }
155   double oneOverC = 1.0/c;                        155   double oneOverC = 1.0/c;
156   pp *= oneOverC;                                 156   pp *= oneOverC;
157   ee *= oneOverC;                                 157   ee *= oneOverC;
158   return *this;                                   158   return *this;
159 } /* w /= c */                                    159 } /* w /= c */
160                                                   160 
161 HepLorentzVector operator / (const HepLorentzV    161 HepLorentzVector operator / (const HepLorentzVector & w, double c) {
162 //  if (c == 0) {                                 162 //  if (c == 0) {
163 //    std::cerr << "HepLorentzVector::operator    163 //    std::cerr << "HepLorentzVector::operator /() - "
164 //      << "Attempt to do LorentzVector / 0 --    164 //      << "Attempt to do LorentzVector / 0 -- \n"
165 //      << "division by zero would produce inf    165 //      << "division by zero would produce infinite or NAN components"
166 //      << std::endl;                             166 //      << std::endl;
167 //  }                                             167 //  }
168   double oneOverC = 1.0/c;                        168   double oneOverC = 1.0/c;
169   return HepLorentzVector (w.getV() * oneOverC    169   return HepLorentzVector (w.getV() * oneOverC,
170                         w.getT() * oneOverC);     170                         w.getT() * oneOverC);
171 } /* LV = w / c */                                171 } /* LV = w / c */
172                                                   172 
173 Hep3Vector HepLorentzVector::boostVector() con    173 Hep3Vector HepLorentzVector::boostVector() const {
174   if (ee == 0) {                                  174   if (ee == 0) {
175     if (pp.mag2() == 0) {                         175     if (pp.mag2() == 0) {
176       return Hep3Vector(0,0,0);                   176       return Hep3Vector(0,0,0);
177     } else {                                      177     } else {
178       std::cerr << "HepLorentzVector::boostVec    178       std::cerr << "HepLorentzVector::boostVector() - "
179         << "boostVector computed for LorentzVe    179         << "boostVector computed for LorentzVector with t=0 -- infinite result"
180         << std::endl;                             180         << std::endl;
181       return pp/ee;                               181       return pp/ee;
182     }                                             182     }
183   }                                               183   }
184   if (restMass2() <= 0) {                         184   if (restMass2() <= 0) {
185     std::cerr << "HepLorentzVector::boostVecto    185     std::cerr << "HepLorentzVector::boostVector() - "
186       << "boostVector computed for a non-timel    186       << "boostVector computed for a non-timelike LorentzVector " << std::endl;
187         // result will make analytic sense but    187         // result will make analytic sense but is physically meaningless
188   }                                               188   }
189   return pp * (1./ee);                            189   return pp * (1./ee);
190 } /* boostVector */                               190 } /* boostVector */
191                                                   191 
192                                                   192 
193 HepLorentzVector & HepLorentzVector::boostX (d    193 HepLorentzVector & HepLorentzVector::boostX (double bbeta){
194   double b2 = bbeta*bbeta;                        194   double b2 = bbeta*bbeta;
195   if (b2 >= 1) {                                  195   if (b2 >= 1) {
196     std::cerr << "HepLorentzVector::boostX() -    196     std::cerr << "HepLorentzVector::boostX() - "
197       << "boost along X with beta >= 1 (speed     197       << "boost along X with beta >= 1 (speed of light) -- \n"
198       << "no boost done" << std::endl;            198       << "no boost done" << std::endl;
199   } else {                                        199   } else {
200     double ggamma = std::sqrt(1./(1-b2));         200     double ggamma = std::sqrt(1./(1-b2));
201     double tt = ee;                               201     double tt = ee;
202     ee = ggamma*(ee + bbeta*pp.getX());           202     ee = ggamma*(ee + bbeta*pp.getX());
203     pp.setX(ggamma*(pp.getX() + bbeta*tt));       203     pp.setX(ggamma*(pp.getX() + bbeta*tt));
204   }                                               204   }
205   return *this;                                   205   return *this;
206 } /* boostX */                                    206 } /* boostX */
207                                                   207 
208 HepLorentzVector & HepLorentzVector::boostY (d    208 HepLorentzVector & HepLorentzVector::boostY (double bbeta){
209   double b2 = bbeta*bbeta;                        209   double b2 = bbeta*bbeta;
210   if (b2 >= 1) {                                  210   if (b2 >= 1) {
211     std::cerr << "HepLorentzVector::boostY() -    211     std::cerr << "HepLorentzVector::boostY() - "
212       << "boost along Y with beta >= 1 (speed     212       << "boost along Y with beta >= 1 (speed of light) -- \n"
213       << "no boost done" << std::endl;            213       << "no boost done" << std::endl;
214   } else {                                        214   } else {
215     double ggamma = std::sqrt(1./(1-b2));         215     double ggamma = std::sqrt(1./(1-b2));
216     double tt = ee;                               216     double tt = ee;
217     ee = ggamma*(ee + bbeta*pp.getY());           217     ee = ggamma*(ee + bbeta*pp.getY());
218     pp.setY(ggamma*(pp.getY() + bbeta*tt));       218     pp.setY(ggamma*(pp.getY() + bbeta*tt));
219   }                                               219   }
220   return *this;                                   220   return *this;
221 } /* boostY */                                    221 } /* boostY */
222                                                   222 
223 HepLorentzVector & HepLorentzVector::boostZ (d    223 HepLorentzVector & HepLorentzVector::boostZ (double bbeta){
224   double b2 = bbeta*bbeta;                        224   double b2 = bbeta*bbeta;
225   if (b2 >= 1) {                                  225   if (b2 >= 1) {
226     std::cerr << "HepLorentzVector::boostZ() -    226     std::cerr << "HepLorentzVector::boostZ() - "
227       << "boost along Z with beta >= 1 (speed     227       << "boost along Z with beta >= 1 (speed of light) -- \n"
228       << "no boost done" << std::endl;            228       << "no boost done" << std::endl;
229   } else {                                        229   } else {
230     double ggamma = std::sqrt(1./(1-b2));         230     double ggamma = std::sqrt(1./(1-b2));
231     double tt = ee;                               231     double tt = ee;
232     ee = ggamma*(ee + bbeta*pp.getZ());           232     ee = ggamma*(ee + bbeta*pp.getZ());
233     pp.setZ(ggamma*(pp.getZ() + bbeta*tt));       233     pp.setZ(ggamma*(pp.getZ() + bbeta*tt));
234   }                                               234   }
235   return *this;                                   235   return *this;
236 } /* boostZ */                                    236 } /* boostZ */
237                                                   237 
238 double HepLorentzVector::setTolerance ( double    238 double HepLorentzVector::setTolerance ( double tol ) {
239 // Set the tolerance for two LorentzVectors to    239 // Set the tolerance for two LorentzVectors to be considered near each other
240   double oldTolerance (tolerance);                240   double oldTolerance (tolerance);
241   tolerance = tol;                                241   tolerance = tol;
242   return oldTolerance;                            242   return oldTolerance;
243 }                                                 243 }
244                                                   244 
245 double HepLorentzVector::getTolerance ( ) {       245 double HepLorentzVector::getTolerance ( ) {
246 // Get the tolerance for two LorentzVectors to    246 // Get the tolerance for two LorentzVectors to be considered near each other
247   return tolerance;                               247   return tolerance;
248 }                                                 248 }
249                                                   249 
250 }  // namespace CLHEP                             250 }  // namespace CLHEP
251                                                   251