Main Page | Class Hierarchy | Compound List | File List | Compound Members | File Members

ParticleFitObject.C

Go to the documentation of this file.
00001 
00018 #include "jbltools/kinfit/ParticleFitObject.h"
00019 #include "cernlib.h"
00020 
00021 #include <iostream>
00022 #include <cassert>
00023 #include <cmath>
00024 using std::isfinite;
00025 
00026 
00027 ParticleFitObject::ParticleFitObject() {
00028   for (int ilocal = 0; ilocal < NPAR; ++ilocal) globalParNum[ilocal] = -1;
00029 }
00030 
00031 ParticleFitObject::~ParticleFitObject()
00032 {}
00033 
00034 bool ParticleFitObject::setParam (int ilocal, double par_, 
00035                                     bool measured_, bool fixed_) {
00036   assert (ilocal >= 0 && ilocal < NPAR);
00037   if (measured[ilocal] != measured_ || fixed[ilocal] != fixed_) invalidateCache();
00038   measured[ilocal] = measured_;
00039   fixed[ilocal] = fixed_;
00040   return setParam (ilocal, par_);
00041 };  
00042 
00043 bool ParticleFitObject::setParam (int ilocal, double par_ ) {
00044   if (!isfinite(par_)) return false;
00045   assert (ilocal >= 0 && ilocal < NPAR);
00046   if (par[ilocal] == par_) return true;
00047   invalidateCache();
00048   par[ilocal] = par_;
00049   return true;
00050 };  
00051 bool ParticleFitObject::setMParam (int ilocal, double mpar_ ) {
00052   if (!isfinite(mpar_)) return false;
00053   assert (ilocal >= 0 && ilocal < NPAR);
00054   if (mpar[ilocal] == mpar_) return true;
00055   invalidateCache();
00056   mpar[ilocal] = mpar_;
00057   return true;
00058 }; 
00059 bool ParticleFitObject::setError (int ilocal, double err_) {
00060   if (!isfinite(err_)) return false;
00061   assert (ilocal >= 0 && ilocal < NPAR);
00062   invalidateCache();
00063   covinvvalid = false;
00064   cov[ilocal][ilocal] = err_*err_;
00065   return true;
00066 }
00067 
00068 bool ParticleFitObject::setCov (int ilocal, int jlocal, double cov_) {
00069   if (!isfinite(cov_)) return false;
00070   assert (ilocal >= 0 && ilocal < NPAR);
00071   assert (jlocal >= 0 && jlocal < NPAR);
00072   invalidateCache();
00073   covinvvalid = false;
00074   cov[ilocal][jlocal] = cov[jlocal][ilocal] = cov_;
00075   return true;
00076 }
00077 bool ParticleFitObject::setMass (double mass_) {
00078   if (!isfinite(mass_)) return false;
00079   if (mass == mass_) return true;
00080   invalidateCache();
00081   mass = std::abs(mass_);
00082   return true;
00083 }
00084 bool ParticleFitObject::fixParam (int ilocal, bool fix) {
00085   assert (ilocal >= 0 && ilocal < NPAR);
00086   return fixed [ilocal] = fix;
00087 }
00088 
00089 bool ParticleFitObject::setGlobalParNum (int ilocal, int iglobal) {
00090   if (ilocal < 0 || ilocal >= NPAR) return false;
00091   globalParNum[ilocal] = iglobal;
00092   return true;
00093 }
00094 int  ParticleFitObject::getGlobalParNum(int ilocal) const {
00095   if (ilocal < 0 || ilocal >= NPAR) return -1;
00096   return globalParNum[ilocal];
00097 }
00098 
00099 double ParticleFitObject::getParam (int ilocal) const {
00100   assert (ilocal >= 0 && ilocal < NPAR);
00101   return par[ilocal];
00102 }
00103 double ParticleFitObject::getMParam (int ilocal) const {
00104   assert (ilocal >= 0 && ilocal < NPAR);
00105   return mpar[ilocal];
00106 }
00107 
00108 double ParticleFitObject::getError (int ilocal) const {
00109   assert (ilocal >= 0 && ilocal < NPAR);
00110   return std::sqrt(cov[ilocal][ilocal]);
00111 }
00112 double ParticleFitObject::getCov (int ilocal, int jlocal) const {
00113   assert (ilocal >= 0 && ilocal < NPAR);
00114   assert (jlocal >= 0 && jlocal < NPAR);
00115   return cov[ilocal][jlocal];
00116 }
00117 bool ParticleFitObject::isParamMeasured (int ilocal) const {
00118   assert (ilocal >= 0 && ilocal < NPAR);
00119   return measured[ilocal];
00120 }
00121 
00122 bool ParticleFitObject::isParamFixed (int ilocal) const {
00123   assert (ilocal >= 0 && ilocal < NPAR);
00124   return fixed[ilocal];
00125 }
00126     
00127 std::ostream&  ParticleFitObject::print4Vector(std::ostream& os) const {
00128   os << "[" << getE() << ", " << getPx() << ", " 
00129       << getPy() << ", "  << getPz() << "]";
00130   return os;
00131 }
00132 
00133 std::ostream&  ParticleFitObject::print (std::ostream& os) const {
00134   printParams(os);
00135   os << " => ";
00136   print4Vector(os);
00137   return os;
00138 }
00139 bool ParticleFitObject::calculateCovInv() const {
00140   int n = getNPar();
00141   int idim = 0;
00142   for (int i = 0; i < n; ++i) {
00143     if (isParamMeasured (i)) {
00144       idim = i;
00145       for (int j = 0; j < n; ++j) {
00146         covinv[i][j] = isParamMeasured (j) ? cov[i][j] : 0;
00147       }
00148     }
00149     else {
00150       for (int j = 0; j < n; ++j) {
00151         covinv[i][j] = static_cast<double>(i == j);
00152       }
00153     }
00154   }
00155   int ierr = (idim == 0) ? 0 : dsinv(idim, &covinv[0][0], NPAR);
00156   if (ierr != 0) {
00157     std::cerr << "ParticleFitObject::calculateCovInv: Error "
00158               << ierr << " from dsinv! Object " << getName() << std::endl;
00159     // printCov (std::cerr);         
00160   }
00161   return covinvvalid = (ierr == 0);
00162 }
00163 
00164 
00165 double ParticleFitObject::getChi2() const {
00166   if (!covinvvalid) calculateCovInv();
00167   if (!covinvvalid) return -1;
00168   double chi2 = 0;
00169   static double resid[NPAR];
00170   static bool chi2contr[NPAR];
00171   for (int i = 0; i < getNPar(); ++i) {
00172     resid[i] = par[i]-mpar[i];
00173     if (chi2contr[i] = isParamMeasured(i) && !isParamFixed(i)) {
00174       chi2 += resid[i]*covinv[i][i]*resid[i];
00175       for (int j = 0; j < i; ++j) {
00176         if (chi2contr[j]) chi2 += 2*(resid[i])*covinv[i][j]*(resid[j]);
00177       }
00178     }
00179   }
00180   return chi2;
00181 }
00182 
00183 double ParticleFitObject::getDChi2DParam(int ilocal) const {
00184   assert (ilocal >= 0 && ilocal < NPAR);
00185   if (isParamFixed(ilocal) || !isParamMeasured(ilocal)) return 0;
00186   if (!covinvvalid) calculateCovInv();
00187   if (!covinvvalid) return 0;
00188   double result = 0;
00189   for (int jlocal = 0; jlocal < getNPar(); jlocal++) 
00190     if (!isParamFixed(jlocal) && isParamMeasured(jlocal))
00191       result += covinv[ilocal][jlocal]*(par[jlocal]-measured[jlocal]);
00192   return 2*result;
00193 }
00194     
00195 double ParticleFitObject::getD2Chi2DParam2(int ilocal, int jlocal) const {
00196   assert (ilocal >= 0 && ilocal < NPAR);
00197   assert (jlocal >= 0 && jlocal < NPAR);
00198   if (isParamFixed(ilocal) || !isParamMeasured(ilocal) && 
00199       isParamFixed(jlocal) || !isParamMeasured(jlocal))
00200     return 0;
00201   if (!covinvvalid) calculateCovInv();
00202   if (!covinvvalid) return 0;
00203   return 2*covinv[ilocal][jlocal];
00204 }
00205       
00206 void ParticleFitObject::addToGlobalChi2DerMatrix (int idim, double *M) const {
00207   if (!covinvvalid) calculateCovInv();
00208   if (!covinvvalid) return;
00209   for (int ilocal = 0; ilocal < getNPar(); ++ilocal) {
00210     if (!isParamFixed(ilocal) && isParamMeasured(ilocal)) {
00211       int iglobal = getGlobalParNum (ilocal);
00212       assert (iglobal >= 0 && iglobal < idim);
00213       int ioffs = idim*iglobal;
00214       for (int jlocal = 0; jlocal < getNPar(); ++jlocal) {
00215         if (!isParamFixed(jlocal) && isParamMeasured(jlocal)) {
00216           int jglobal = getGlobalParNum (jlocal);
00217           assert (jglobal >= 0 && jglobal < idim);
00218           M[ioffs+jglobal] += 2*covinv[ilocal][jlocal];
00219         }
00220       }
00221     }
00222   }
00223 }
00224 
00225 
00226 void ParticleFitObject::addToGlobCov(double *globCov, int idim) const {
00227   for (int ilocal = 0; ilocal < getNPar(); ++ilocal) {
00228     if (!isParamFixed(ilocal) && isParamMeasured(ilocal)) {
00229       int iglobal = getGlobalParNum (ilocal);
00230       assert (iglobal >= 0 && iglobal < idim);
00231       int ioffs = idim*iglobal;
00232       for (int jlocal = 0; jlocal < getNPar(); ++jlocal) {
00233         if (!isParamFixed(jlocal) && isParamMeasured(jlocal)) {
00234           int jglobal = getGlobalParNum (jlocal);
00235           assert (jglobal >= 0 && jglobal < idim);
00236           globCov[ioffs+jglobal] += cov[ilocal][jlocal];
00237         }
00238       }
00239     }
00240   }
00241 }
00242 
00243 

Generated on Fri Sep 14 17:38:21 2007 for Kinfit by doxygen 1.3.2