/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Mdc/DedxCorrecSvc/DedxCorrecSvc-00-02-52/src/Dedx_Helix.cxx

Go to the documentation of this file.

//
// $Id: Dedx_Helix.cxx,v 1.3 2011/02/23 00:36:21 maqm Exp $
//
// $Log: Dedx_Helix.cxx,v $
// Revision 1.3  2011/02/23 00:36:21  maqm
// slc5-x64-gcc43
//
// Revision 1.2  2011/02/17 11:09:39  maqm
// x64-slc5-gcc43
//
// Revision 1.1  2008/06/30 01:43:15  xcao
// change Helix into Dedx_Helix
//
// Revision 1.2  2008/04/18 03:36:05  max
// *** empty log message ***
//
// Revision 1.1.1.1  2006/01/06 08:25:05  wangdy
// first import of this package
//
// Revision 1.2  2005/04/20 13:56:14  wangdy
// codes updated for Bes detector ,see recent changeLog
//
// Revision 1.1.1.1  2005/03/17 13:31:56  wangdy
// first import this package
//
// Revision 1.19  2002/01/03 11:05:06  katayama
// Point3D and other header files are cleaned
//
// Revision 1.18  2001/05/07 21:06:37  yiwasaki
// <float.h> included for linux
//
// Revision 1.17  2001/04/25 02:55:39  yiwasaki
// cache m_ac[5] added
//
// Revision 1.16  2001/04/18 12:09:19  katayama
// optimized
//
// Revision 1.15  2000/01/26 10:22:53  yiwasaki
// copy operator bug fix(reported by M.Yokoyama)
//
// Revision 1.14  1999/11/23 10:29:22  yiwasaki
// static cosnt double Helix::ConstantAlpha added
//
// Revision 1.13  1999/06/15 01:49:41  yiwasaki
// constructor bug fix
//
// Revision 1.12  1999/06/15 01:44:53  yiwasaki
// minor change
//
// Revision 1.11  1999/06/06 07:34:27  yiwasaki
// i/o functions to set/get mag. field added
//
// Revision 1.10  1999/05/11 23:26:59  yiwasaki
// option to ignore error calculations added
//
// Revision 1.9  1998/11/06 08:51:19  katayama
// protect 0div
//
// Revision 1.8  1998/07/08 04:35:50  jtanaka
// add some members to the constructor by Iwasaki-san and Ozaki-san
//
// Revision 1.7  1998/06/18 10:27:20  katayama
// Added several new functions from Tanaka san
//
// Revision 1.6  1998/02/24 01:07:59  yiwasaki
// minor fix
//
// Revision 1.5  1998/02/24 00:31:52  yiwasaki
// bug fix
//
// Revision 1.4  1998/02/20 02:22:23  yiwasaki
// for new helix class
//
// Revision 1.3  1998/01/22 08:38:07  hitoshi
// Fixed a bug in fid cal. (found at Osaka).
//
// Revision 1.2  1997/09/19 00:35:52  katayama
// Added Id and Log
//
//
//
//   Class Helix 
//
//   Author      Date         comments
//   Y.Ohnishi   03/01/1997   original version
//   Y.Ohnishi   06/03/1997   updated
//   Y.Iwasaki   17/02/1998   BFILED removed, func. name changed, func. added
//   J.Tanaka    06/12/1998   add some utilities.
//   Y.Iwasaki   07/07/1998   cache added to speed up
//

#include <math.h>
#include <float.h>
#include <iostream>
#include "DedxCorrecSvc/Dedx_Helix.h"
#include "CLHEP/Matrix/Matrix.h"
#include "GaudiKernel/StatusCode.h"
#include "GaudiKernel/IInterface.h"
#include "GaudiKernel/Kernel.h"
#include "GaudiKernel/Service.h"
#include "GaudiKernel/Kernel.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/SvcFactory.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/IMessageSvc.h"
#include "GaudiKernel/IPartPropSvc.h"

#include "GaudiKernel/SmartDataPtr.h"


const double
M_PI2 = 2. * M_PI;

const double
M_PI4 = 4. * M_PI;

const double
M_PI8 = 8. * M_PI;

const double
Dedx_Helix::ConstantAlpha = -333.564095;

Dedx_Helix::Dedx_Helix(const HepPoint3D & pivot,
             const HepVector & a,
             const HepSymMatrix & Ea)
: //m_bField(10.0),
  //m_alpha(333.564095),
  m_pivot(pivot),
  m_a(a),
  m_matrixValid(true),
  m_Ea(Ea) {
    //ISvcLocator * SvcLocator = Gaudi::SvcLocator(); 
    StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF); 
    if(scmgn!=StatusCode::SUCCESS) { 
     // log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq; 
     std::cout<< "Unable to open Magnetic field service"<<std::endl;
    }
    m_bField = 10000*(m_pmgnIMF->getReferField());
    m_alpha = 10000. / 2.99792458 / m_bField;
    // m_alpha = 222.376063;
    updateCache();
}

Dedx_Helix::Dedx_Helix(const HepPoint3D & pivot,
             const HepVector & a)
: //m_bField(10.0),
  //m_alpha(333.564095),
  m_pivot(pivot),
  m_a(a),
  m_matrixValid(false),
  m_Ea(HepSymMatrix(5,0)) {
    StatusCode scmgn = Gaudi::svcLocator()->service ("MagneticFieldSvc",m_pmgnIMF);
    if(scmgn!=StatusCode::SUCCESS) {
     // log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
     std::cout<< "Unable to open Magnetic field service"<<std::endl;
    }
    m_bField = 10000*(m_pmgnIMF->getReferField());
    m_alpha = 10000. / 2.99792458 / m_bField;
    //cout<<"m_bField is = "<<m_bField<<"        m_alpha is = "<<m_alpha<<endl;
    // m_alpha = 222.376063;
    updateCache();
}

Dedx_Helix::Dedx_Helix(const HepPoint3D & position,
             const Hep3Vector & momentum,
             double charge)
: //m_bField(10.0),
  //m_alpha(333.564095),
  m_pivot(position),
  m_a(HepVector(5,0)),
  m_matrixValid(false),
  m_Ea(HepSymMatrix(5,0)) {
    StatusCode scmgn = Gaudi::svcLocator()->service ("MagneticFieldSvc",m_pmgnIMF);
    if(scmgn!=StatusCode::SUCCESS) {
     // log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
     std::cout<< "Unable to open Magnetic field service"<<std::endl;
    }
    m_bField = 10000*(m_pmgnIMF->getReferField());
    m_alpha = 10000. / 2.99792458 / m_bField;
    
    m_a[0] = 0.;
    m_a[1] = fmod(atan2(- momentum.x(), momentum.y())
                  + M_PI4, M_PI2);
    m_a[3] = 0.;
    double perp(momentum.perp());
    if (perp != 0.0) {
        m_a[2] = charge / perp;
        m_a[4] = momentum.z() / perp; 
    }
    else {
        m_a[2] = charge * (DBL_MAX);
        if (momentum.z() >= 0) {
            m_a[4] = (DBL_MAX);
        } else {
            m_a[4] = -(DBL_MAX);
        }
    }
    // m_alpha = 222.376063;
    updateCache();
}

Dedx_Helix::~Dedx_Helix() {
}

HepPoint3D
Dedx_Helix::x(double phi) const {
    //
    // Calculate position (x,y,z) along helix.
    //   
    // x = x0 + dr * cos(phi0) + (alpha / kappa) * (cos(phi0) - cos(phi0+phi))
    // y = y0 + dr * sin(phi0) + (alpha / kappa) * (sin(phi0) - sin(phi0+phi))
    // z = z0 + dz             - (alpha / kappa) * tan(lambda) * phi
    //
   // m_ac[0] = m_ac[0] *(-1);
    // cout<<"m_bField is = "<<m_bField<<"        m_alpha is = "<<m_alpha<<endl;
  // cout<<"m_pivot is = "<<m_pivot<<"       m_ac is = ("<<m_ac[0]<<" , "<<m_ac[1]<<" , "<<m_ac[2]<<" , "<<m_ac[3]<<" , "<<m_ac[4]<<" )"<<endl;
  //cout<<"m_cp is = "<<m_cp<<"     "<<"   m_sp is =  "<<m_sp<<"      m_r is =  "<<m_r<<endl;  
    //cout<<" m_r + m_ac[0] = "<<m_ac[0] +m_r<<endl;
    double x = m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] +phi));
    double y = m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] +phi));
    double z = m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi;

    return HepPoint3D(x, y, z);
}

double *
Dedx_Helix::x(double phi, double p[3]) const {
    //
    // Calculate position (x,y,z) along helix.
    //   
    // x = x0 + dr * cos(phi0) + (alpha / kappa) * (cos(phi0) - cos(phi0+phi))
    // y = y0 + dr * sin(phi0) + (alpha / kappa) * (sin(phi0) - sin(phi0+phi))
    // z = z0 + dz             - (alpha / kappa) * tan(lambda) * phi
    //

    p[0] = m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] + phi));
    p[1] = m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] + phi));
    p[2] = m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi;

    return p;
}

HepPoint3D
Dedx_Helix::x(double phi, HepSymMatrix & Ex) const {
    double x = m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] +phi));
    double y = m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] +phi));
    double z = m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi;

    //
    //   Calculate position error matrix.
    //   Ex(phi) = (@x/@a)(Ea)(@x/@a)^T, phi is deflection angle to specify the
    //   point to be calcualted.
    //
    // HepMatrix dXDA(3, 5, 0);
    // dXDA = delXDelA(phi);
    // Ex.assign(dXDA * m_Ea * dXDA.T());

    if (m_matrixValid) Ex = m_Ea.similarity(delXDelA(phi));
    else               Ex = m_Ea;

    return HepPoint3D(x, y, z);
}

Hep3Vector
Dedx_Helix::momentum(double phi) const {
    // 
    // Calculate momentum.
    //
    // Pt = | 1/kappa | (GeV/c)
    //
    // Px = -Pt * sin(phi0 + phi) 
    // Py =  Pt * cos(phi0 + phi)
    // Pz =  Pt * tan(lambda)
    //

    double pt = fabs(m_pt);
    double px = - pt * sin(m_ac[1] + phi);
    double py =   pt * cos(m_ac[1] + phi);
    double pz =   pt * m_ac[4];

    return Hep3Vector(px, py, pz);
}

Hep3Vector
Dedx_Helix::momentum(double phi, HepSymMatrix & Em) const {
    // 
    // Calculate momentum.
    //
    // Pt = | 1/kappa | (GeV/c)
    //
    // Px = -Pt * sin(phi0 + phi) 
    // Py =  Pt * cos(phi0 + phi)
    // Pz =  Pt * tan(lambda)
    //

    double pt = fabs(m_pt);
    double px = - pt * sin(m_ac[1] + phi);
    double py =   pt * cos(m_ac[1] + phi);
    double pz =   pt * m_ac[4];

    if (m_matrixValid) Em = m_Ea.similarity(delMDelA(phi));
    else               Em = m_Ea;

    return Hep3Vector(px, py, pz);
}

HepLorentzVector 
Dedx_Helix::momentum(double phi, double mass) const {
    // 
    // Calculate momentum.
    //
    // Pt = | 1/kappa | (GeV/c)
    //
    // Px = -Pt * sin(phi0 + phi) 
    // Py =  Pt * cos(phi0 + phi)
    // Pz =  Pt * tan(lambda)
    //
    // E  = sqrt( 1/kappa/kappa * (1+tan(lambda)*tan(lambda)) + mass*mass )

    double pt = fabs(m_pt);
    double px = - pt * sin(m_ac[1] + phi);
    double py =   pt * cos(m_ac[1] + phi);
    double pz =   pt * m_ac[4];
    double E  =   sqrt(pt*pt*(1.+m_ac[4]*m_ac[4])+mass*mass);

    return HepLorentzVector(px, py, pz, E);
}


HepLorentzVector 
Dedx_Helix::momentum(double phi, double mass, HepSymMatrix & Em) const {
    // 
    // Calculate momentum.
    //
    // Pt = | 1/kappa | (GeV/c)
    //
    // Px = -Pt * sin(phi0 + phi) 
    // Py =  Pt * cos(phi0 + phi)
    // Pz =  Pt * tan(lambda)
    //
    // E  = sqrt( 1/kappa/kappa * (1+tan(lambda)*tan(lambda)) + mass*mass )

    double pt = fabs(m_pt);
    double px = - pt * sin(m_ac[1] + phi);
    double py =   pt * cos(m_ac[1] + phi);
    double pz =   pt * m_ac[4];
    double E  =   sqrt(pt*pt*(1.+m_ac[4]*m_ac[4])+mass*mass);

    if (m_matrixValid) Em = m_Ea.similarity(del4MDelA(phi,mass));
    else               Em = m_Ea;

    return HepLorentzVector(px, py, pz, E);
}

HepLorentzVector 
Dedx_Helix::momentum(double phi,
                double mass,
                HepPoint3D & x,
                HepSymMatrix & Emx) const {
  // 
  // Calculate momentum.
  //
  // Pt = | 1/kappa | (GeV/c)
  //
  // Px = -Pt * sin(phi0 + phi) 
  // Py =  Pt * cos(phi0 + phi)
  // Pz =  Pt * tan(lambda)
  //
  // E  = sqrt( 1/kappa/kappa * (1+tan(lambda)*tan(lambda)) + mass*mass )
  
  double pt = fabs(m_pt);
  double px = - pt * sin(m_ac[1] + phi);
  double py =   pt * cos(m_ac[1] + phi);
  double pz =   pt * m_ac[4];
  double E  = sqrt(pt * pt * (1. + m_ac[4] * m_ac[4]) + mass * mass);

  x.setX(m_pivot.x() + m_ac[0] * m_cp + m_r * (m_cp - cos(m_ac[1] + phi)));
  x.setY(m_pivot.y() + m_ac[0] * m_sp + m_r * (m_sp - sin(m_ac[1] + phi)));
  x.setZ(m_pivot.z() + m_ac[3] - m_r * m_ac[4] * phi);

  if (m_matrixValid) Emx = m_Ea.similarity(del4MXDelA(phi,mass));
  else               Emx = m_Ea;
  
  return HepLorentzVector(px, py, pz, E);
}


const HepPoint3D &
Dedx_Helix::pivot(const HepPoint3D & newPivot) {
    const double & dr    = m_ac[0];
    const double & phi0  = m_ac[1];
    const double & kappa = m_ac[2];
    const double & dz    = m_ac[3];
    const double & tanl  = m_ac[4];

    double rdr = dr + m_r;
    double phi = fmod(phi0 + M_PI4, M_PI2);
    double csf0 = cos(phi);
    double snf0 = (1. - csf0) * (1. + csf0);
    snf0 = sqrt((snf0 > 0.) ? snf0 : 0.);
    if(phi > M_PI) snf0 = - snf0;

    double xc = m_pivot.x() + rdr * csf0;
    double yc = m_pivot.y() + rdr * snf0;
    double csf, snf;
    if(m_r != 0.0) {
      csf = (xc - newPivot.x()) / m_r;
      snf = (yc - newPivot.y()) / m_r;
      double anrm = sqrt(csf * csf + snf * snf);
      if(anrm != 0.0) {
        csf /= anrm;
        snf /= anrm;
        phi = atan2(snf, csf);
      } else {
        csf = 1.0;
        snf = 0.0;
        phi = 0.0;
      }
    } else {
      csf = 1.0;
      snf = 0.0;
      phi = 0.0;
    }
    double phid = fmod(phi - phi0 + M_PI8, M_PI2);
    if(phid > M_PI) phid = phid - M_PI2;
    double drp = (m_pivot.x() + dr * csf0 + m_r * (csf0 - csf) - newPivot.x())
        * csf
        + (m_pivot.y() + dr * snf0 + m_r * (snf0 - snf) - newPivot.y()) * snf;
    double dzp = m_pivot.z() + dz - m_r * tanl * phid - newPivot.z();

    HepVector ap(5);
    ap[0] = drp;
    ap[1] = fmod(phi + M_PI4, M_PI2);
    ap[2] = kappa;
    ap[3] = dzp;
    ap[4] = tanl;

    //    if (m_matrixValid) m_Ea.assign(delApDelA(ap) * m_Ea * delApDelA(ap).T());
    if (m_matrixValid) m_Ea = m_Ea.similarity(delApDelA(ap));

    m_a = ap;
    m_pivot = newPivot;

    //...Are these needed?...iw...
    updateCache();
    return m_pivot;
}

void
Dedx_Helix::set(const HepPoint3D & pivot,
           const HepVector & a,
           const HepSymMatrix & Ea) {
    m_pivot = pivot;
    m_a = a;
    m_Ea = Ea;
    m_matrixValid = true;
    updateCache();
}

Dedx_Helix &
Dedx_Helix::operator = (const Dedx_Helix & i) {
    if (this == & i) return * this;

    m_bField = i.m_bField;
    m_alpha = i.m_alpha;
    m_pivot = i.m_pivot;
    m_a = i.m_a;
    m_Ea = i.m_Ea;
    m_matrixValid = i.m_matrixValid;

    m_center = i.m_center;
    m_cp = i.m_cp;
    m_sp = i.m_sp;
    m_pt = i.m_pt;
    m_r  = i.m_r;
    m_ac[0] = i.m_ac[0];
    m_ac[1] = i.m_ac[1];
    m_ac[2] = i.m_ac[2];
    m_ac[3] = i.m_ac[3];
    m_ac[4] = i.m_ac[4];

    return * this;
}

void
Dedx_Helix::updateCache(void) {
    //
    //   Calculate Helix center( xc, yc ).
    //   
    //   xc = x0 + (dr + (alpha / kappa)) * cos(phi0)  (cm)
    //   yc = y0 + (dr + (alpha / kappa)) * sin(phi0)  (cm)
    //

    m_ac[0] = m_a[0];
    m_ac[1] = m_a[1];
    m_ac[2] = m_a[2];
    m_ac[3] = m_a[3];
    m_ac[4] = m_a[4];

    m_cp = cos(m_ac[1]);
    m_sp = sin(m_ac[1]);
    if (m_ac[2] != 0.0) {
        m_pt = 1. / m_ac[2];
        m_r = m_alpha / m_ac[2];
    }
    else {
        m_pt = (DBL_MAX);
        m_r = (DBL_MAX);
    }

    double x = m_pivot.x() + (m_ac[0] + m_r) * m_cp;
    double y = m_pivot.y() + (m_ac[0] + m_r) * m_sp;
    m_center.setX(x);
    m_center.setY(y);
    m_center.setZ(0.);
}

HepMatrix
Dedx_Helix::delApDelA(const HepVector & ap) const {
    //
    //   Calculate Jacobian (@ap/@a)
    //   Vector ap is new helix parameters and a is old helix parameters. 
    //

    HepMatrix dApDA(5,5,0);

    const double & dr    = m_ac[0];
    const double & phi0  = m_ac[1];
    const double & cpa   = m_ac[2];
    const double & dz    = m_ac[3];
    const double & tnl   = m_ac[4];

    double drp   = ap[0];
    double phi0p = ap[1];
    double cpap  = ap[2];
    double dzp   = ap[3];
    double tnlp  = ap[4];

    double rdr   = m_r + dr;
    double rdrpr;
    if ((m_r + drp) != 0.0) {
      rdrpr = 1. / (m_r + drp);
    } else {
      rdrpr = (DBL_MAX);
    }
    // double csfd  = cos(phi0)*cos(phi0p) + sin(phi0)*sin(phi0p); 
    // double snfd  = cos(phi0)*sin(phi0p) - sin(phi0)*cos(phi0p); 
    double csfd  = cos(phi0p - phi0);
    double snfd  = sin(phi0p - phi0);
    double phid  = fmod(phi0p - phi0 + M_PI8, M_PI2);
    if (phid > M_PI) phid = phid - M_PI2;

    dApDA[0][0]  =  csfd;
    dApDA[0][1]  =  rdr*snfd;
    if(cpa!=0.0) {
      dApDA[0][2]  =  (m_r/cpa)*( 1.0 - csfd );
    } else {
      dApDA[0][2]  = (DBL_MAX);
    }

    dApDA[1][0]  = - rdrpr*snfd;
    dApDA[1][1]  = rdr*rdrpr*csfd;
    if(cpa!=0.0) {
      dApDA[1][2]  = (m_r/cpa)*rdrpr*snfd;
    } else {
      dApDA[1][2]  = (DBL_MAX);
    }
    
    dApDA[2][2]  = 1.0;

    dApDA[3][0]  = m_r*rdrpr*tnl*snfd;
    dApDA[3][1]  = m_r*tnl*(1.0 - rdr*rdrpr*csfd);
    if(cpa!=0.0) {
      dApDA[3][2]  = (m_r/cpa)*tnl*(phid - m_r*rdrpr*snfd);
    } else {
      dApDA[3][2]  = (DBL_MAX);
    }
    dApDA[3][3]  = 1.0;
    dApDA[3][4]  = - m_r*phid;

    dApDA[4][4] = 1.0;
 
    return dApDA;
}

HepMatrix
Dedx_Helix::delXDelA(double phi) const {
    //
    //   Calculate Jacobian (@x/@a)
    //   Vector a is helix parameters and phi is internal parameter
    //   which specifys the point to be calculated for Ex(phi).
    //

    HepMatrix dXDA(3,5,0);

    const double & dr      = m_ac[0];
    const double & phi0    = m_ac[1];
    const double & cpa     = m_ac[2];
    const double & dz      = m_ac[3];
    const double & tnl     = m_ac[4];

    double cosf0phi = cos(phi0 + phi);
    double sinf0phi = sin(phi0 + phi);

    dXDA[0][0]     = m_cp;
    dXDA[0][1]     = - dr * m_sp + m_r * (- m_sp + sinf0phi); 
    if(cpa!=0.0) {
      dXDA[0][2]     = - (m_r / cpa) * (m_cp - cosf0phi);
    } else {
      dXDA[0][2] = (DBL_MAX);
    }
    // dXDA[0][3]     = 0.0;
    // dXDA[0][4]     = 0.0;
 
    dXDA[1][0]     = m_sp;
    dXDA[1][1]     = dr * m_cp + m_r * (m_cp - cosf0phi);
    if(cpa!=0.0) {
      dXDA[1][2]     = - (m_r / cpa) * (m_sp - sinf0phi);
    } else {
      dXDA[1][2] = (DBL_MAX);
    }
    // dXDA[1][3]     = 0.0;
    // dXDA[1][4]     = 0.0;

    // dXDA[2][0]     = 0.0;
    // dXDA[2][1]     = 0.0;
    if(cpa!=0.0) {
      dXDA[2][2]     = (m_r / cpa) * tnl * phi;
    } else {
      dXDA[2][2] = (DBL_MAX);
    }
    dXDA[2][3]     = 1.0;
    dXDA[2][4]     = - m_r * phi;

    return dXDA;
}



HepMatrix
Dedx_Helix::delMDelA(double phi) const {
    //
    //   Calculate Jacobian (@m/@a)
    //   Vector a is helix parameters and phi is internal parameter.
    //   Vector m is momentum.
    //

    HepMatrix dMDA(3,5,0);

    const double & phi0 = m_ac[1];
    const double & cpa  = m_ac[2];
    const double & tnl  = m_ac[4];

    double cosf0phi = cos(phi0+phi);
    double sinf0phi = sin(phi0+phi);

    double rho;
    if(cpa != 0.)rho = 1./cpa;
    else rho = (DBL_MAX);

    double charge = 1.;
    if(cpa < 0.)charge = -1.;

    dMDA[0][1] = -fabs(rho)*cosf0phi;
    dMDA[0][2] = charge*rho*rho*sinf0phi;
 
    dMDA[1][1] = -fabs(rho)*sinf0phi;
    dMDA[1][2] = -charge*rho*rho*cosf0phi;

    dMDA[2][2] = -charge*rho*rho*tnl;
    dMDA[2][4] = fabs(rho);

    return dMDA;
}


HepMatrix
Dedx_Helix::del4MDelA(double phi, double mass) const {
    //
    //   Calculate Jacobian (@4m/@a)
    //   Vector a  is helix parameters and phi is internal parameter.
    //   Vector 4m is 4 momentum.
    //

    HepMatrix d4MDA(4,5,0);

    double phi0 = m_ac[1];
    double cpa  = m_ac[2];
    double tnl  = m_ac[4];

    double cosf0phi = cos(phi0+phi);
    double sinf0phi = sin(phi0+phi);

    double rho;
    if(cpa != 0.)rho = 1./cpa;
    else rho = (DBL_MAX);

    double charge = 1.;
    if(cpa < 0.)charge = -1.;

    double E = sqrt(rho*rho*(1.+tnl*tnl)+mass*mass);

    d4MDA[0][1] = -fabs(rho)*cosf0phi;
    d4MDA[0][2] = charge*rho*rho*sinf0phi;
 
    d4MDA[1][1] = -fabs(rho)*sinf0phi;
    d4MDA[1][2] = -charge*rho*rho*cosf0phi;

    d4MDA[2][2] = -charge*rho*rho*tnl;
    d4MDA[2][4] = fabs(rho);

    if (cpa != 0.0 && E != 0.0) {
      d4MDA[3][2] = (-1.-tnl*tnl)/(cpa*cpa*cpa*E);
      d4MDA[3][4] = tnl/(cpa*cpa*E);
    } else {
      d4MDA[3][2] = (DBL_MAX);
      d4MDA[3][4] = (DBL_MAX);
    }
    return d4MDA;
}


HepMatrix
Dedx_Helix::del4MXDelA(double phi, double mass) const {
    //
    //   Calculate Jacobian (@4mx/@a)
    //   Vector a  is helix parameters and phi is internal parameter.
    //   Vector 4xm is 4 momentum and position.
    //

    HepMatrix d4MXDA(7,5,0);

    const double & dr      = m_ac[0];
    const double & phi0    = m_ac[1];
    const double & cpa     = m_ac[2];
    const double & dz      = m_ac[3];
    const double & tnl     = m_ac[4];

    double cosf0phi = cos(phi0+phi);
    double sinf0phi = sin(phi0+phi);

    double rho;
    if(cpa != 0.)rho = 1./cpa;
    else rho = (DBL_MAX);

    double charge = 1.;
    if(cpa < 0.)charge = -1.;

    double E = sqrt(rho * rho * (1. + tnl * tnl) + mass * mass);

    d4MXDA[0][1] = - fabs(rho) * cosf0phi;
    d4MXDA[0][2] = charge * rho * rho * sinf0phi;
 
    d4MXDA[1][1] = - fabs(rho) * sinf0phi;
    d4MXDA[1][2] = - charge * rho * rho * cosf0phi;

    d4MXDA[2][2] = - charge * rho * rho * tnl;
    d4MXDA[2][4] = fabs(rho);

    if (cpa != 0.0 && E != 0.0) {
      d4MXDA[3][2] = (- 1. - tnl * tnl) / (cpa * cpa * cpa * E);
      d4MXDA[3][4] = tnl / (cpa * cpa * E);
    } else {
      d4MXDA[3][2] = (DBL_MAX);
      d4MXDA[3][4] = (DBL_MAX);
    }
    
    d4MXDA[4][0] = m_cp;
    d4MXDA[4][1] = - dr * m_sp + m_r * (- m_sp + sinf0phi);
    if (cpa != 0.0) {
      d4MXDA[4][2] = - (m_r / cpa) * (m_cp - cosf0phi);
    } else {
      d4MXDA[4][2] = (DBL_MAX);
    }
     
    d4MXDA[5][0] = m_sp;
    d4MXDA[5][1] = dr * m_cp + m_r * (m_cp - cosf0phi);
    if (cpa != 0.0) {
      d4MXDA[5][2] = - (m_r / cpa) * (m_sp - sinf0phi);

      d4MXDA[6][2] = (m_r / cpa) * tnl * phi;
    } else {
      d4MXDA[5][2] = (DBL_MAX);

      d4MXDA[6][2] = (DBL_MAX);
    }

    d4MXDA[6][3] = 1.;
    d4MXDA[6][4] = - m_r * phi;

    return d4MXDA;
}

void
Dedx_Helix::ignoreErrorMatrix(void) {
    m_matrixValid = false;
    m_Ea *= 0.;
}

---