StableChargedParticle.C

Go to the documentation of this file.

#include"StableChargedParticle.h"
#include"cernlib.h"
#include <cmath>
#include <cassert>

using std::sqrt;
using std::sin;
using std::cos;

// constructor
StableChargedParticle::StableChargedParticle(double kappa, double phi0, double theta, 
                                             double dca, double z0,
                                             double bfield_, double mass_)
  : bfield (bfield_), mass (mass_)
                                             {
  setParam (0, kappa, true);
  setParam (1, phi0, true);
  setParam (2, theta, true);
  setParam (3, dca, true);
  setParam (4, z0, true);
  setParam (5, 0, false);   // start tracklength
  setMass (m);
  for (int ilocal = 0; ilocal < NPAR; ++ilocal) globalParNum[ilocal] = -1;
  invalidateCache();
}

// destructor
StableChargedParticle::~StableChargedParticle() {}

// these depend on actual parametrisation!


    
void   StableChargedParticle::addToDerivatives (double der[],
                                       int idim,
                                       double pxfact, 
                                       double pyfact, 
                                       double pzfact, 
                                       double efact) const {
  assert (0);                                     
  int i_E     = globalParNum[0];
  int i_theta = globalParNum[1];
  int i_phi   = globalParNum[2];
  assert (i_E     >= 0 && i_E     < idim);
  assert (i_theta >= 0 && i_theta < idim);
  assert (i_phi   >= 0 && i_phi   < idim);
  
  if (!cachevalid) updateCache();
  // for numerical accuracy, add up derivatives first,
  // then add them to global vector
  double der_E = 0;
  double der_theta = 0;
  double der_phi = 0;
  
  if (pxfact != 0) {
    der_E     += pxfact*dpxdE;
    der_theta += pxfact*dpydtheta;
    der_phi   -= pxfact*py;
  }
  if (pyfact != 0) {
    der_E     += pyfact*dpydE;
    der_theta += pyfact*dpydtheta;
    der_phi   += pyfact*px;
  }
  if (pzfact != 0) {
    der_E     += pzfact*dpdE*ctheta;
    der_theta -= pzfact*pt;
  }
  der_E     += efact;
  
  der[i_E]     += der_E;
  der[i_theta] += der_theta;
  der[i_phi]   += der_phi;
}
    
void   StableChargedParticle::addTo2ndDerivatives (double der2[],
                                          int idim,
                                          double pxfact, 
                                          double pyfact, 
                                          double pzfact, 
                                          double efact) const {
  int i_E  = globalParNum[0];
  int i_th = globalParNum[1];
  int i_ph= globalParNum[2];
  assert (i_E  >= 0 && i_E  < idim);
  assert (i_th >= 0 && i_th < idim);
  assert (i_ph >= 0 && i_ph < idim);
  
  if (!cachevalid) updateCache();
  // for numerical accuracy, add up derivatives first,
  // then add them to global vector
  double der_EE   = 0;
  double der_Eth  = 0;
  double der_Eph  = 0;
  double der_thth = 0;
  double der_thph = 0;
  double der_phph = 0;
  
  double d2pdE2 = (mass != 0) ? -mass*mass/(p*p*p) : 0;
  double d2ptdE2 = d2pdE2*stheta;
  
  if (pxfact != 0) {
    der_EE   += pxfact*d2ptdE2*cphi;
    der_Eth  += pxfact*dpzdE*cphi;
    der_Eph  -= pxfact*dpydE;
    der_thth -= pxfact*px;
    der_thph -= pxfact*dpydtheta;
    der_phph -= pxfact*px;
  }
  if (pyfact != 0) {
    der_EE   += pyfact*d2ptdE2*sphi;
    der_Eth  += pyfact*dpzdE*sphi;
    der_Eph  += pyfact*dpxdE;
    der_thth -= pyfact*py;
    der_thph += pyfact*dpxdtheta;
    der_phph -= pyfact*py;
  }
  if (pzfact != 0) {
    der_EE   += pzfact*d2pdE2*ctheta;
    der_Eth  -= pzfact*dptdE;
    der_thth -= pzfact*pz;
  }
  
  der2[idim*i_E+i_E]   += der_EE;
  der2[idim*i_E+i_th]  += der_Eth;
  der2[idim*i_E+i_ph]  += der_Eph;
  der2[idim*i_th+i_E]  += der_Eth;
  der2[idim*i_th+i_th] += der_thth;
  der2[idim*i_th+i_ph] += der_thph;
  der2[idim*i_ph+i_E]  += der_Eph;
  der2[idim*i_ph+i_th] += der_thph;
  der2[idim*i_ph+i_ph] += der_phph;
}


    
double StableChargedParticle::getChi2() const {
  if (!cachevalid) updateCache();
  return chi2;
}
      
void StableChargedParticle::addToGlobalChi2DerMatrix (int idim, double *M) const {
  assert (0);
  // Only diagonal terms d^2 chi^2 / dx_i dx_i are nonzero
  for (int ilocal = 0; ilocal < getNPar(); ++ilocal) {
    int iglobal = getGlobalParNum (ilocal);
    assert (iglobal >= 0 && iglobal < idim);
    M[(idim+1)*iglobal] += 2./std::pow (err[ilocal], 2);
  }
}
  
      
void StableChargedParticle::addToGlobalDerMatrix (int idim, double c, double *M) const {
  assert (0);
  // add second derivatives to global matrix
  if (c == 0) return;
  
  
}
  
double StableChargedParticle::getDChi2DParam(int ilocal) const {
  assert (ilocal >= 0 && ilocal < NPAR);
  if (ilocal >= NMEAS) return 0;
  if (!cachevalid) updateCache();
  double result = 0;
  for (int i = 0; i<NMEAS; i++) result += 2*covinv[ilocal][i]*resid[i];
  return ;
}
    
double StableChargedParticle::getD2Chi2DParam2(int ilocal1, int ilocal2) const {
  assert (ilocal1 >= 0 && ilocal1 < NPAR);
  assert (ilocal2 >= 0 && ilocal2 < NPAR);
  if (ilocal1 >= NMEAS || ilocal2 >= NMEAS ) return 0;
  return 2*covinv[ilocal1][ilocal2];
}

void StableChargedParticle::setCovMatrix(double values[NMEAS*(NMEAS+1)/2]) {
  for (int i = 0; i < NMEAS; ++i) {
    for (int j = 0; j < NMEAS; ++j) {
      covmatrix[i][j] = 0.;
    }
  } 
  int k=0;   
  for (int i = 0; i < NMEAS; ++i) {
    covinv[i][i] = covmatrix[i][i] = std::pow(values[k++], 2);
  }
  for (int i = 0; i < 4; ++i) {
    for (int j = i+1; j < NMEAS; ++j) {
      covinv[i][j] = covinv[j][i] = covmatrix[i][j] = covmatrix[j][i] = values[k++];
    }
  }
  int ifail = dsinv (NMEAS, covinv, NMEAS);
  assert (ifail == 0); 
}

void StableChargedParticle::addToGlobCov(double *cov, int idim) const {
  assert(0);
  for (int ilocal = 0; ilocal < NPAR; ilocal++) {
    int iglobal = getGlobalParNum(ilocal);
    if (iglobal >= 0) cov[iglobal*(1+idim)] += err[ilocal]*err[ilocal];
//     std::cout << "cov[" << iglobal*(1+idim) << "] = " << cov[iglobal*(1+idim)]
//               << " for idim = " << idim << " and ilocal = " << ilocal << std::endl;
  }
}

void StableChargedParticle::invalidateCache() const {
  cachevalid = false;
}

void StableChargedParticle::updateCache() const {
  double& kappa = par[0];
  double& phi0  = par[1];
  double& theta = par[2];
  double& dca   = par[3];
  double& z0    = par[4];
  double& s1    = par[5];
  
  sphi0 = std::sin(phi0);
  cphi0 = std::cos(phi0);
  r = (kappa != 0) ? 1/kappa : 0;
  dcamir = dca - r;
  cottheta = 1/std::tan(theta);
  cos2theta = std::pow (std::cos(theta), 2);
  
  setVertex (v1, dv1, 5, s1);
  
  chi2 = 0;
  for (int i = 0; i < NMEAS; ++i) {
    resid[i] = par[i]-measured[i];
    chi2 += resid[i]*covinv[i][i]*resid[i];
    for (int j = 0; j < i; ++j) {
      chi2 += 2*(resid[i])*covinv[i][j]*(resid[j])
    }
  }
  cachevalid = true;
}
void StableChargedParticle::setVertex (FitThreeVector& v, FitThreeVector dv[NPAR], 
                                       int n, double s) const {
  double& kappa = par[0];
  double& dca   = par[3];
  double& z0    = par[4];
  
  if (std::abs (kappa*s) < 1.e-6) {
    double ssq = s*s;
    double kssq = 0.5*kappa*ssq;
    v.setValues (sphi0*( dca + kssq) - cphi0*s,
                 cphi0*(-dca + kssq) + sphi0*s,
                 z0 + s*cottheta);
    dv[0].setValues (-0.5*sphi0*ssq,
                      0.5*cphi0*ssq, 
                     0);          
  }
  else {
    double kappas = kappa*s;
    double phi0plkappas = phi0 + kappas;
    double si = std::sin(phi0plkappas);
    double co = std::cos(phi0plkappas);
    v.setValues ( dcamir*sphi0 + r*si,
                 -dcamir*cphi0 - r*co,
                 z0 + s*cottheta);
    dv[0].setValues (co*kappas + sphi0 -si,
                     si*kappas - cphi0 +co, 
                     0);          
  }
  dv[1].setValues (-v.getY(),
                   -v.getX(), 
                   0);          
  dv[2].setValues (0, 0, - s/cos2theta);          
  dv[3].setValues (sphi0, -cphi0, 0);          
  dv[4].setValues (0, 0, 1);  
  assert (n<NPAR);        
  dv[n].setValues (co, si, cottheta);          
}

---