Math/GenVector/eta.h

0001 // @(#)root/mathcore:$Id$
0002 // Authors: W. Brown, M. Fischler, L. Moneta    2005
0003
0004  /**********************************************************************
0005   *                                                                    *
0006   * Copyright (c) 2005 , FNAL MathLib Team                             *
0007   *                                                                    *
0008   *                                                                    *
0009   **********************************************************************/
0010
0011
0012 // Header source file for function calculating eta
0013 //
0014 // Created by: Lorenzo Moneta  at 14 Jun 2007
0015
0016
0017 #ifndef ROOT_Math_GenVector_eta
0018 #define ROOT_Math_GenVector_eta  1
0019
0020 #include "Math/GenVector/etaMax.h"
0021
0022
0023 #include <limits>
0024 #include <cmath>
0025
0026
0027 namespace ROOT {
0028
0029   namespace Math {
0030
0031      namespace Impl {
0032
0033     /**
0034         Calculate eta given rho and zeta.
0035         This formula is faster than the standard calculation (below) from log(tan(theta/2)
0036         but one has to be careful when rho is much smaller than z (large eta values)
0037         Formula is  eta = log( zs + sqrt(zs^2 + 1) )  where zs = z/rho
0038
0039         For large value of z_scaled (tan(theta) ) one can approximate the sqrt via a Taylor expansion
0040         We do the approximation of the sqrt if the numerical error is of the same order of second term of
0041         the sqrt.expansion:
0042         eps > 1/zs^4   =>   zs > 1/(eps^0.25)
0043
0044         When rho == 0 we use etaMax (see definition in etaMax.h)
0045
0046      */
0047         template<typename Scalar>
0048         inline Scalar Eta_FromRhoZ(Scalar rho, Scalar z) {
0049            if (rho > 0) {
0050
0051               // value to control Taylor expansion of sqrt
0052               static const Scalar big_z_scaled = pow(std::numeric_limits<Scalar>::epsilon(), static_cast<Scalar>(-.25));
0053
0054               Scalar z_scaled = z/rho;
0055               using std::fabs;
0056               if (std::fabs(z_scaled) < big_z_scaled) {
0057                  using std::sqrt;
0058                  using std::log;
0059                  return log(z_scaled + std::sqrt(z_scaled * z_scaled + 1.0));
0060               } else {
0061                  // apply correction using first order Taylor expansion of sqrt
0062                  using std::log;
0063                  return z > 0 ? log(2.0 * z_scaled + 0.5 / z_scaled) : -log(-2.0 * z_scaled);
0064               }
0065            }
0066            // case vector has rho = 0
0067            else if (z==0) {
0068               return 0;
0069            }
0070            else if (z>0) {
0071               return z + etaMax<Scalar>();
0072            }
0073            else {
0074               return z - etaMax<Scalar>();
0075            }
0076
0077         }
0078
0079
0080         /**
0081            Implementation of eta from -log(tan(theta/2)).
0082            This is convenient when theta is already known (for example in a polar coorindate system)
0083         */
0084         template<typename Scalar>
0085         inline Scalar Eta_FromTheta(Scalar theta, Scalar r) {
0086            Scalar tanThetaOver2 = tan(theta / 2.);
0087            if (tanThetaOver2 == 0) {
0088               return r + etaMax<Scalar>();
0089            }
0090            else if (tanThetaOver2 > std::numeric_limits<Scalar>::max()) {
0091               return -r - etaMax<Scalar>();
0092            }
0093            else {
0094               using std::log;
0095               return -log(tanThetaOver2);
0096            }
0097
0098         }
0099
0100      } // end namespace Impl
0101
0102   } // namespace Math
0103
0104 } // namespace ROOT
0105
0106
0107 #endif /* ROOT_Math_GenVector_etaMax  */