Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/externals/clhep/src/LorentzVectorC.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/LorentzVectorC.cc (Version 11.3.0) and /externals/clhep/src/LorentzVectorC.cc (Version 11.0.p4)


  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 the HepLorent      6 // This is the implementation of the HepLorentzVector class:
  7 // Those methods originating with ZOOM dealing      7 // Those methods originating with ZOOM dealing with comparison (other than
  8 // isSpaceLike, isLightlike, isTimelike, which      8 // isSpaceLike, isLightlike, isTimelike, which are in the main part.)
  9 //                                                  9 //
 10 // 11/29/05 mf in deltaR, replaced the direct      10 // 11/29/05 mf in deltaR, replaced the direct subtraction 
 11 // pp.phi() - w.getV().phi() with pp.deltaRPhi     11 // pp.phi() - w.getV().phi() with pp.deltaRPhi(w.getV()) which behaves 
 12 // correctly across the 2pi boundary.              12 // correctly across the 2pi boundary.
 13                                                    13 
 14 #include "CLHEP/Vector/LorentzVector.h"            14 #include "CLHEP/Vector/LorentzVector.h"
 15                                                    15 
 16 #include <cmath>                                   16 #include <cmath>
 17 #include <iostream>                                17 #include <iostream>
 18                                                    18 
 19 namespace CLHEP  {                                 19 namespace CLHEP  {
 20                                                    20 
 21 //-***********                                     21 //-***********
 22 // Comparisons                                     22 // Comparisons
 23 //-***********                                     23 //-***********
 24                                                    24 
 25 int HepLorentzVector::compare (const HepLorent     25 int HepLorentzVector::compare (const HepLorentzVector & w) const {
 26   if       ( ee > w.ee ) {                         26   if       ( ee > w.ee ) {
 27     return 1;                                      27     return 1;
 28   } else if ( ee < w.ee ) {                        28   } else if ( ee < w.ee ) {
 29     return -1;                                     29     return -1;
 30   } else {                                         30   } else {
 31     return ( pp.compare(w.pp) );                   31     return ( pp.compare(w.pp) );
 32   }                                                32   }
 33 } /* Compare */                                    33 } /* Compare */
 34                                                    34 
 35 bool HepLorentzVector::operator > (const HepLo     35 bool HepLorentzVector::operator > (const HepLorentzVector & w) const {
 36         return (compare(w)  > 0);                  36         return (compare(w)  > 0);
 37 }                                                  37 }
 38 bool HepLorentzVector::operator < (const HepLo     38 bool HepLorentzVector::operator < (const HepLorentzVector & w) const {
 39         return (compare(w)  < 0);                  39         return (compare(w)  < 0);
 40 }                                                  40 }
 41 bool HepLorentzVector::operator>= (const HepLo     41 bool HepLorentzVector::operator>= (const HepLorentzVector & w) const {
 42         return (compare(w) >= 0);                  42         return (compare(w) >= 0);
 43 }                                                  43 }
 44 bool HepLorentzVector::operator<= (const HepLo     44 bool HepLorentzVector::operator<= (const HepLorentzVector & w) const {
 45         return (compare(w) <= 0);                  45         return (compare(w) <= 0);
 46 }                                                  46 }
 47                                                    47 
 48 //-********                                        48 //-********
 49 // isNear                                          49 // isNear
 50 // howNear                                         50 // howNear
 51 //-********                                        51 //-********
 52                                                    52 
 53 bool HepLorentzVector::isNear(const HepLorentz     53 bool HepLorentzVector::isNear(const HepLorentzVector & w, 
 54             double epsilon) const {                54             double epsilon) const {
 55   double limit = std::fabs(pp.dot(w.pp));          55   double limit = std::fabs(pp.dot(w.pp));
 56   limit += .25*((ee+w.ee)*(ee+w.ee));              56   limit += .25*((ee+w.ee)*(ee+w.ee));
 57   limit *= epsilon*epsilon;                        57   limit *= epsilon*epsilon;
 58   double delta = (pp - w.pp).mag2();               58   double delta = (pp - w.pp).mag2();
 59   delta +=  (ee-w.ee)*(ee-w.ee);                   59   delta +=  (ee-w.ee)*(ee-w.ee);
 60   return (delta <= limit );                        60   return (delta <= limit );
 61 } /* isNear() */                                   61 } /* isNear() */
 62                                                    62 
 63 double HepLorentzVector::howNear(const HepLore     63 double HepLorentzVector::howNear(const HepLorentzVector & w) const {
 64   double wdw = std::fabs(pp.dot(w.pp)) + .25*(     64   double wdw = std::fabs(pp.dot(w.pp)) + .25*((ee+w.ee)*(ee+w.ee));
 65   double delta = (pp - w.pp).mag2() + (ee-w.ee     65   double delta = (pp - w.pp).mag2() + (ee-w.ee)*(ee-w.ee);
 66   if ( (wdw > 0) && (delta < wdw)  ) {             66   if ( (wdw > 0) && (delta < wdw)  ) {
 67     return std::sqrt (delta/wdw);                  67     return std::sqrt (delta/wdw);
 68   } else if ( (wdw == 0) && (delta == 0) ) {       68   } else if ( (wdw == 0) && (delta == 0) ) {
 69     return 0;                                      69     return 0;
 70   } else {                                         70   } else {
 71     return 1;                                      71     return 1;
 72   }                                                72   }
 73 } /* howNear() */                                  73 } /* howNear() */
 74                                                    74 
 75 //-*********                                       75 //-*********
 76 // isNearCM                                        76 // isNearCM
 77 // howNearCM                                       77 // howNearCM
 78 //-*********                                       78 //-*********
 79                                                    79 
 80 bool HepLorentzVector::isNearCM                    80 bool HepLorentzVector::isNearCM
 81       (const HepLorentzVector & w, double epsi     81       (const HepLorentzVector & w, double epsilon) const {
 82                                                    82 
 83   double tTotal = (ee + w.ee);                     83   double tTotal = (ee + w.ee);
 84   Hep3Vector vTotal (pp + w.pp);                   84   Hep3Vector vTotal (pp + w.pp);
 85   double vTotal2 = vTotal.mag2();                  85   double vTotal2 = vTotal.mag2();
 86                                                    86 
 87   if ( vTotal2 >= tTotal*tTotal ) {                87   if ( vTotal2 >= tTotal*tTotal ) {
 88     // Either one or both vectors are spacelik     88     // Either one or both vectors are spacelike, or the dominant T components
 89     // are in opposite directions.  So boostin     89     // are in opposite directions.  So boosting and testing makes no sense;
 90     // but we do consider two exactly equal ve     90     // but we do consider two exactly equal vectors to be equal in any frame,
 91     // even if they are spacelike and can't be     91     // even if they are spacelike and can't be boosted to a CM frame.
 92     return (*this == w);                           92     return (*this == w);
 93   }                                                93   }
 94                                                    94 
 95   if ( vTotal2 == 0 ) {  // no boost needed!       95   if ( vTotal2 == 0 ) {  // no boost needed!
 96     return (isNear(w, epsilon));                   96     return (isNear(w, epsilon));
 97   }                                                97   }
 98                                                    98 
 99   // Find the boost to the CM frame.  We know      99   // Find the boost to the CM frame.  We know that the total vector is timelike.
100                                                   100 
101   double tRecip = 1./tTotal;                      101   double tRecip = 1./tTotal;
102   Hep3Vector bboost ( vTotal * (-tRecip) );       102   Hep3Vector bboost ( vTotal * (-tRecip) );
103                                                   103 
104         //-| Note that you could do pp/t and n    104         //-| Note that you could do pp/t and not be terribly inefficient since
105         //-| SpaceVector/t itself takes 1/t an    105         //-| SpaceVector/t itself takes 1/t and multiplies.  The code here saves
106         //-| a redundant check for t=0.           106         //-| a redundant check for t=0.
107                                                   107 
108   // Boost both vectors.  Since we have the sa    108   // Boost both vectors.  Since we have the same boost, there is no need
109   // to repeat the beta and gamma calculation;    109   // to repeat the beta and gamma calculation; and there is no question
110   // about beta >= 1.  That is why we don't ju    110   // about beta >= 1.  That is why we don't just call w.boosted().
111                                                   111 
112   double b2 = vTotal2*tRecip*tRecip;              112   double b2 = vTotal2*tRecip*tRecip;
113                                                   113 
114   double ggamma = std::sqrt(1./(1.-b2));          114   double ggamma = std::sqrt(1./(1.-b2));
115   double boostDotV1 = bboost.dot(pp);             115   double boostDotV1 = bboost.dot(pp);
116   double gm1_b2 = (ggamma-1)/b2;                  116   double gm1_b2 = (ggamma-1)/b2;
117                                                   117 
118   HepLorentzVector w1 ( pp   + ((gm1_b2)*boost    118   HepLorentzVector w1 ( pp   + ((gm1_b2)*boostDotV1+ggamma*ee) * bboost,
119                      ggamma * (ee + boostDotV1    119                      ggamma * (ee + boostDotV1) );
120                                                   120 
121   double boostDotV2 = bboost.dot(w.pp);           121   double boostDotV2 = bboost.dot(w.pp);
122   HepLorentzVector w2 ( w.pp + ((gm1_b2)*boost    122   HepLorentzVector w2 ( w.pp + ((gm1_b2)*boostDotV2+ggamma*w.ee) * bboost,
123                      ggamma * (w.ee + boostDot    123                      ggamma * (w.ee + boostDotV2) );
124                                                   124 
125   return (w1.isNear(w2, epsilon));                125   return (w1.isNear(w2, epsilon));
126                                                   126 
127 } /* isNearCM() */                                127 } /* isNearCM() */
128                                                   128 
129 double HepLorentzVector::howNearCM(const HepLo    129 double HepLorentzVector::howNearCM(const HepLorentzVector & w) const {
130                                                   130 
131   double tTotal = (ee + w.ee);                    131   double tTotal = (ee + w.ee);
132   Hep3Vector vTotal (pp + w.pp);                  132   Hep3Vector vTotal (pp + w.pp);
133   double vTotal2 = vTotal.mag2();                 133   double vTotal2 = vTotal.mag2();
134                                                   134 
135   if ( vTotal2 >= tTotal*tTotal ) {               135   if ( vTotal2 >= tTotal*tTotal ) {
136     // Either one or both vectors are spacelik    136     // Either one or both vectors are spacelike, or the dominant T components
137     // are in opposite directions.  So boostin    137     // are in opposite directions.  So boosting and testing makes no sense;
138     // but we do consider two exactly equal ve    138     // but we do consider two exactly equal vectors to be equal in any frame,
139     // even if they are spacelike and can't be    139     // even if they are spacelike and can't be boosted to a CM frame.
140     if (*this == w) {                             140     if (*this == w) {
141       return 0;                                   141       return 0;
142     } else {                                      142     } else {
143       return 1;                                   143       return 1;
144     }                                             144     }
145   }                                               145   }
146                                                   146 
147   if ( vTotal2 == 0 ) {  // no boost needed!      147   if ( vTotal2 == 0 ) {  // no boost needed!
148     return (howNear(w));                          148     return (howNear(w));
149   }                                               149   }
150                                                   150 
151   // Find the boost to the CM frame.  We know     151   // Find the boost to the CM frame.  We know that the total vector is timelike.
152                                                   152 
153   double tRecip = 1./tTotal;                      153   double tRecip = 1./tTotal;
154   Hep3Vector bboost ( vTotal * (-tRecip) );       154   Hep3Vector bboost ( vTotal * (-tRecip) );
155                                                   155 
156         //-| Note that you could do pp/t and n    156         //-| Note that you could do pp/t and not be terribly inefficient since
157         //-| SpaceVector/t itself takes 1/t an    157         //-| SpaceVector/t itself takes 1/t and multiplies.  The code here saves
158         //-| a redundant check for t=0.           158         //-| a redundant check for t=0.
159                                                   159 
160   // Boost both vectors.  Since we have the sa    160   // Boost both vectors.  Since we have the same boost, there is no need
161   // to repeat the beta and gamma calculation;    161   // to repeat the beta and gamma calculation; and there is no question
162   // about beta >= 1.  That is why we don't ju    162   // about beta >= 1.  That is why we don't just call w.boosted().
163                                                   163 
164   double b2 = vTotal2*tRecip*tRecip;              164   double b2 = vTotal2*tRecip*tRecip;
165 //  if ( b2 >= 1 ) {      // NaN-proofing         165 //  if ( b2 >= 1 ) {      // NaN-proofing
166 //    std::cerr << "HepLorentzVector::howNearC    166 //    std::cerr << "HepLorentzVector::howNearCM() - "
167 //  << "boost vector in howNearCM appears to b    167 //  << "boost vector in howNearCM appears to be tachyonic" << std::endl;
168 //  }                                             168 //  }
169   double ggamma = std::sqrt(1./(1.-b2));          169   double ggamma = std::sqrt(1./(1.-b2));
170   double boostDotV1 = bboost.dot(pp);             170   double boostDotV1 = bboost.dot(pp);
171   double gm1_b2 = (ggamma-1)/b2;                  171   double gm1_b2 = (ggamma-1)/b2;
172                                                   172 
173   HepLorentzVector w1 ( pp   + ((gm1_b2)*boost    173   HepLorentzVector w1 ( pp   + ((gm1_b2)*boostDotV1+ggamma*ee) * bboost,
174                      ggamma * (ee + boostDotV1    174                      ggamma * (ee + boostDotV1) );
175                                                   175 
176   double boostDotV2 = bboost.dot(w.pp);           176   double boostDotV2 = bboost.dot(w.pp);
177   HepLorentzVector w2 ( w.pp + ((gm1_b2)*boost    177   HepLorentzVector w2 ( w.pp + ((gm1_b2)*boostDotV2+ggamma*w.ee) * bboost,
178                      ggamma * (w.ee + boostDot    178                      ggamma * (w.ee + boostDotV2) );
179                                                   179 
180   return (w1.howNear(w2));                        180   return (w1.howNear(w2));
181                                                   181 
182 } /* howNearCM() */                               182 } /* howNearCM() */
183                                                   183 
184 //-************                                   184 //-************
185 // deltaR                                         185 // deltaR
186 // isParallel                                     186 // isParallel
187 // howParallel                                    187 // howParallel
188 // howLightlike                                   188 // howLightlike
189 //-************                                   189 //-************
190                                                   190 
191 double HepLorentzVector::deltaR ( const HepLor    191 double HepLorentzVector::deltaR ( const HepLorentzVector & w ) const {
192                                                   192 
193   double a = eta() - w.eta();                     193   double a = eta() - w.eta();
194   double b = pp.deltaPhi(w.getV());               194   double b = pp.deltaPhi(w.getV());
195                                                   195 
196   return std::sqrt ( a*a + b*b );                 196   return std::sqrt ( a*a + b*b );
197                                                   197 
198 } /* deltaR */                                    198 } /* deltaR */
199                                                   199 
200 // If the difference (in the Euclidean norm) o    200 // If the difference (in the Euclidean norm) of the normalized (in Euclidean
201 // norm) 4-vectors is small, then those 4-vect    201 // norm) 4-vectors is small, then those 4-vectors are considered nearly
202 // parallel.                                      202 // parallel.
203                                                   203 
204 bool HepLorentzVector::isParallel (const HepLo    204 bool HepLorentzVector::isParallel (const HepLorentzVector & w, double epsilon) const {
205   double norm = euclideanNorm();                  205   double norm = euclideanNorm();
206   double wnorm = w.euclideanNorm();               206   double wnorm = w.euclideanNorm();
207   if ( norm == 0 ) {                              207   if ( norm == 0 ) {
208     if ( wnorm == 0 ) {                           208     if ( wnorm == 0 ) {
209       return true;                                209       return true;
210     } else {                                      210     } else {
211       return false;                               211       return false;
212     }                                             212     }
213   }                                               213   }
214   if ( wnorm == 0 ) {                             214   if ( wnorm == 0 ) {
215     return false;                                 215     return false;
216   }                                               216   }
217   HepLorentzVector w1 = *this / norm;             217   HepLorentzVector w1 = *this / norm;
218   HepLorentzVector w2 = w / wnorm;                218   HepLorentzVector w2 = w / wnorm;
219   return ( (w1-w2).euclideanNorm2() <= epsilon    219   return ( (w1-w2).euclideanNorm2() <= epsilon*epsilon );
220 } /* isParallel */                                220 } /* isParallel */
221                                                   221 
222                                                   222 
223 double HepLorentzVector::howParallel (const He    223 double HepLorentzVector::howParallel (const HepLorentzVector & w) const {
224                                                   224 
225   double norm = euclideanNorm();                  225   double norm = euclideanNorm();
226   double wnorm = w.euclideanNorm();               226   double wnorm = w.euclideanNorm();
227   if ( norm == 0 ) {                              227   if ( norm == 0 ) {
228     if ( wnorm == 0 ) {                           228     if ( wnorm == 0 ) {
229       return 0;                                   229       return 0;
230     } else {                                      230     } else {
231       return 1;                                   231       return 1;
232     }                                             232     }
233   }                                               233   }
234   if ( wnorm == 0 ) {                             234   if ( wnorm == 0 ) {
235     return 1;                                     235     return 1;
236   }                                               236   }
237                                                   237 
238   HepLorentzVector w1 = *this / norm;             238   HepLorentzVector w1 = *this / norm;
239   HepLorentzVector w2 = w / wnorm;                239   HepLorentzVector w2 = w / wnorm;
240   double x1 = (w1-w2).euclideanNorm();            240   double x1 = (w1-w2).euclideanNorm();
241   return (x1 < 1) ? x1 : 1;                       241   return (x1 < 1) ? x1 : 1;
242                                                   242 
243 } /* howParallel */                               243 } /* howParallel */
244                                                   244 
245 double HepLorentzVector::howLightlike() const     245 double HepLorentzVector::howLightlike() const {
246   double m1 = std::fabs(restMass2());             246   double m1 = std::fabs(restMass2());
247   double twoT2 = 2*ee*ee;                         247   double twoT2 = 2*ee*ee;
248   if (m1 < twoT2) {                               248   if (m1 < twoT2) {
249     return m1/twoT2;                              249     return m1/twoT2;
250   } else {                                        250   } else {
251     return 1;                                     251     return 1;
252   }                                               252   }
253 } /* HowLightlike */                              253 } /* HowLightlike */
254                                                   254 
255 }  // namespace CLHEP                             255 }  // namespace CLHEP
256                                                   256