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

PConstraint.C

Go to the documentation of this file.
00001 
00002 // Class PConstraint
00003 //
00004 // Author: Jenny Boehme, Bennmo List
00005 // Last update: $Date: 2005/01/12 10:11:45 $
00006 //          by: $Author: blist $
00007 // 
00008 // Description: constraint 
00009 // a*sum(px)+b*sum(py)+c*sum(pz)+d*sum(E)=e
00010 //               
00012 
00013 #include "jbltools/kinfit/PConstraint.h"
00014 #include "jbltools/kinfit/ParticleFitObject.h"
00015 
00016 #include<iostream>
00017 #include<cassert>
00018 
00019 using std::cout;
00020 using std::endl;
00021 
00022 PConstraint::PConstraint (double pxfact_, double pyfact_, double pzfact_,
00023                           double efact_, double value_) 
00024     : pxfact (pxfact_),
00025       pyfact (pyfact_),
00026       pzfact (pzfact_),
00027       efact (efact_),
00028       value (value_),
00029       cachevalid(false)
00030 {}
00031 
00032 // destructor
00033 PConstraint::~PConstraint () {}
00034 
00035 // calculate current value of constraint function
00036 double PConstraint::getValue() const {
00037   double totpx = 0;
00038   double totpy = 0;
00039   double totpz = 0;
00040   double totE = 0;
00041   for (unsigned int i = 0; i < fitobjects.size(); i++) {
00042     if (pxfact != 0) totpx += fitobjects[i]->getPx(); 
00043     if (pyfact != 0) totpy += fitobjects[i]->getPy(); 
00044     if (pzfact != 0) totpz += fitobjects[i]->getPz(); 
00045     if (efact  != 0) totE  += fitobjects[i]->getE(); 
00046   }
00047   return pxfact*totpx + pyfact*totpy + pzfact*totpz + efact*totE - value;
00048 }
00049 
00050 // calculate vector/array of derivatives of this contraint 
00051 // w.r.t. to ALL parameters of all fitobjects
00052 // here: d sum(px) /d par(i,j) 
00053 //                      = d sum(px) /d px(i) * d px(i) /d par(i, j)
00054 //                                      =  1 * d px(i) /d par(i, j)
00055 void PConstraint::getDerivatives(int idim, double der[]) const {
00056   for (unsigned int i = 0; i < fitobjects.size(); i++) {
00057     for (int ilocal = 0; ilocal < fitobjects[i]->getNPar(); ilocal++) {
00058       if (!fitobjects[i]->isParamFixed(ilocal)) {
00059         int iglobal = fitobjects[i]->getGlobalParNum (ilocal);
00060         assert (iglobal >= 0 && iglobal < idim);
00061         double d = 0;
00062         if (pxfact != 0) d += pxfact*fitobjects[i]->getDPx (ilocal);
00063         if (pyfact != 0) d += pyfact*fitobjects[i]->getDPy (ilocal);
00064         if (pzfact != 0) d += pzfact*fitobjects[i]->getDPz (ilocal);
00065         if (efact  != 0) d +=  efact*fitobjects[i]->getDE  (ilocal);
00066         der[iglobal] = d;
00067       }
00068     }
00069   }
00070 }
00071   
00072 void PConstraint::add1stDerivativesToMatrix(int idim, double *M) const {
00073 
00074   assert (0);
00075 
00076   assert (M);
00077   int kglobal = getGlobalNum();
00078   assert (kglobal >= 0 && kglobal < idim);
00079   
00080   for (ConstFitObjectIterator i = fitobjects.begin(); i != fitobjects.end(); ++i) {
00081     const ParticleFitObject *fo = *i;
00082     assert (fo);
00083     for (int ilocal = 0; ilocal < fo->getNPar(); ++ilocal) {
00084       if (!fo->isParamFixed(ilocal)) {
00085         int iglobal = fo->getGlobalParNum (ilocal);
00086         assert (iglobal >= 0 && iglobal < idim);
00087         double d = fo->getDPx(ilocal);
00088         M[idim*iglobal+kglobal] += d;
00089         M[idim*kglobal+iglobal] += d;
00090       }
00091     }
00092   }
00093 }
00094   
00095 void PConstraint::add2ndDerivativesToMatrix(int idim, double *M, double lambda) const {
00096 
00097   assert (0);
00098 }
00099     
00100 void PConstraint::addToGlobalDerMatrix (double lambda, int idim, double *M) const {
00101 
00102   assert (0);
00103   // Add lambda*d^2 g / d x_i dx_j to global matrix
00104   
00105   if (lambda == 0) return;
00106   
00107   // d^2 g / (dx_i dx_j) = 
00108   //   = sum_k,l d^2 g/(dpx_k dpx_l) * dpx_k/dx_i dpx_l/dx_j
00109   //     + sum_k dg/dpx_k * d^2 px_k/(dx_i dx_j)
00110   //   = sum_k,l      1              * dpx_k/dx_i dpx_l/dx_j
00111   
00112   // assume here that different 4-vectors always depend on 
00113   // different parameters!
00114   
00115   if (!cachevalid) updateCache();
00116   
00117   int *globalParNum = new int[nparams];
00118   double *der = new double[nparams];
00119   
00120 //   ipar = 0;
00121 //   for (int i = 0; i < fitobjects.size(); i++) {
00122 //     for (int ilocal = 0; ilocal < fitobjects[i]->getNPar(); ilocal++) {
00123 //       int iglobal = fitobjects[i]->getGlobalParNum (ilocal);
00124 //       if (iglobal >= 0) {
00125 //         assert (ipar < nparams);
00126 //         globalParNum[ipar] = iglobal;
00127 //         der[ipar] = fitobjects[i]->getDPx (ilocal);
00128 //         ipar++;
00129 //       }
00130 //   }
00131   
00132   for (int ipar = 0; ipar < nparams; ipar++) {
00133     int iglobal = globalParNum[ipar];
00134     double der_i = der[ipar];
00135     for (int jpar = ipar; jpar < nparams; jpar++) {
00136       int jglobal = globalParNum[ipar];
00137       double der_j = der[jpar];
00138       double l_der_ij = lambda*der_i*der_j;
00139       M[idim*iglobal+jglobal] += l_der_ij;
00140       if (ipar != jpar) M[idim*jglobal+iglobal] += l_der_ij;
00141     }
00142   }       
00143 }
00144 
00145 void PConstraint::invalidateCache() const {
00146   cachevalid = false;
00147 }
00148 
00149 void PConstraint::updateCache() const {
00150   nparams = 0;
00151   for (unsigned int i = 0; i < fitobjects.size(); i++) {
00152     for (int ilocal = 0; ilocal < fitobjects[i]->getNPar(); ilocal++) {
00153       int iglobal = fitobjects[i]->getGlobalParNum (ilocal);
00154       if (!fitobjects[i]->isParamFixed(ilocal)) {
00155         assert (iglobal >= 0);
00156         nparams++;
00157       }
00158     }
00159   }
00160   cachevalid = true;
00161 }
00162   

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