DifNumber.cxx

Go to the documentation of this file.

//--------------------------------------------------------------------------
// File and Version Information:
//  $Id: DifNumber.cxx,v 1.3 2010/03/25 09:55:57 zhangy Exp $
//
// Description:
//  Class Implementation for |DifNumber| 
//      What do i do ?
// Environment:
//  Software developed for the BaBar Detector at the SLAC B-Factory.
//
// Author List:
//  A. Snyder
//
// Copyright Information:
//  Copyright (C) 1996  SLAC
// 
// History:
//  Migration for BESIII MDC
//
//------------------------------------------------------------------------
 
#include "MdcRecoUtil/DifNumber.h"
#include "MdcRecoUtil/DifIndepPar.h"
using std::endl;
//CHANGE using std::ostream;
using std::cout;
 
extern const DifNumber zero(0.0);
extern const DifNumber one(1.0);
 
 
double DifNumber::error(const HepSymMatrix& e)const {
  return sqrt(correlation(*this,e));
}
 
double DifNumber::error()const{
  if(indepPar()==0) {return 0.0;}
  return error(indepPar()->covariance());
}
 
void DifNumber::fetchDerivatives(HepVector& v)const {
  assert(v.num_row()==nPar());
  for(int i=1; i<=nPar(); i++) {v(i)=derivative(i);}
}
 
HepVector DifNumber::derivatives()const{
    HepVector temp(nPar());
    fetchDerivatives(temp);
    return temp;
}
 
double 
DifNumber::correlation(const DifNumber& b,const HepSymMatrix &e) const {
  assert(e.num_col()==nPar());
  assert(e.num_row()==b.nPar());
  double error = 0.;
  for(int i=1; i<=nPar(); i++) {
    for(int j=1; j<=b.nPar(); j++) {
      error+=derivative(i)*e(i,j)*b.derivative(j);
    }
  }
  return error;
}
 
double
DifNumber::correlation(const DifNumber& b)const {
  if(indepPar()==0) return 0.0;
  if(b.indepPar()!=indepPar()) return 0.0;
  return correlation(b,indepPar()->covariance());
}
 
double correlation(const DifNumber& a,const DifNumber& b) // correlation from default matrix
{ 
  return (a.indepPar()==0||b.indepPar()==0||a.indepPar()!=b.indepPar())?0:a.correlation(b,a.indepPar()->covariance()); 
}
 
// FIXME: This function should be inlined, but that would require checking out additional packages...
double correlation(const DifNumber& a,const DifNumber& b,const HepSymMatrix& e) // correlation for specified error
{ return a.correlation(b,e); }
 
void DifNumber::print(/*ostream& o*/)const {
  cout << "number:" << number() << endl;
  cout << "npar:" << nPar() << endl;
  for(int i=1; i<=nPar(); i++) {
    cout << "derivative(" << i << "):" << derivative(i) << endl;
  }
}
 
 
extern DifNumber solveQuad
(const DifNumber& a,        // quadratic term
 const DifNumber& b,        // linear term
 const DifNumber& c,        // const term
 int pref,          // solution preference
 Code& code)            // error code
{
  DifNumber descr=b*b-4.0*a*c;
  if(descr<0.0) {       // solution not real
    code.setFail(1341);
    return DifNumber(0.0);
  }
  if(a.number()==0.0){
    if(b.number()==0.0) {
      code.setFail(1342);
      return DifNumber(0.0);
    }
    code.setSuccess(40);
    return -c/b+a*c/pow(b,3);
  }
  code.setSuccess(40);
  descr=sqrt(descr);
  DifNumber s=-b;
 
  if(pref==+1) {        // positive solution
    s+=descr;
  }else if(pref==-1){       // negative solution
    s-=descr;
  }else if(pref==0) {       // smallest solution
    if(s>0.0) {s-=descr;}else {s+=descr;}
  }else {           // illegal prefrence
    code.setFail(1343);
    return DifNumber(0.0);
  }
  s/=2.0*a;
  return s;
}