KinFitter/src/TAbsFitParticle.cc

0001 // Classname: TAbsFitParticle
0002 // Author: Jan E. Sundermann, Verena Klose (TU Dresden)
0003
0004 //________________________________________________________________
0005 //
0006 // TAbsFitParticle::
0007 // --------------------
0008 //
0009 // Abstract base class for particles to be used with kinematic fitter
0010 //
0011
0012 #include "PhysicsTools/KinFitter/interface/TAbsFitParticle.h"
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014 #include <iostream>
0015 #include <iomanip>
0016 #include "TClass.h"
0017
0018 TAbsFitParticle::TAbsFitParticle()
0019     : TNamed("NoName", "NoTitle"),
0020       _nPar(0),
0021       _u1(),
0022       _u2(),
0023       _u3(),
0024       _covMatrix(),
0025       _covMatrixFit(),
0026       _covMatrixDeltaY(),
0027       _pull(),
0028       _iniparameters(1, 1),
0029       _parameters(1, 1),
0030       _pini(),
0031       _pcurr() {}
0032
0033 TAbsFitParticle::TAbsFitParticle(const TString& name, const TString& title)
0034     : TNamed(name, title),
0035       _nPar(0),
0036       _u1(),
0037       _u2(),
0038       _u3(),
0039       _covMatrix(),
0040       _covMatrixFit(),
0041       _covMatrixDeltaY(),
0042       _pull(),
0043       _iniparameters(1, 1),
0044       _parameters(1, 1),
0045       _pini(),
0046       _pcurr() {}
0047
0048 TAbsFitParticle::~TAbsFitParticle() {}
0049
0050 TString TAbsFitParticle::getInfoString() {
0051   // Collect information to be used for printout
0052
0053   std::stringstream info;
0054   info << std::scientific << std::setprecision(6);
0055
0056   info << "__________________________" << std::endl << std::endl;
0057
0058   info << "OBJ: " << IsA()->GetName() << "\t" << GetName() << "\t" << GetTitle() << std::endl;
0059
0060   info << std::setw(22) << "initial parameters:" << std::setw(5) << " " << std::setw(20)
0061        << "current parameters:" << std::endl;
0062   for (int i = 0; i < _nPar; i++) {
0063     info << "par[" << i << "] = " << std::setw(18) << (*getParIni())(i, 0) << std::setw(20) << (*getParCurr())(i, 0)
0064          << std::endl;
0065   }
0066
0067   info << std::setw(22) << "initial 4vector:" << std::setw(5) << " " << std::setw(20)
0068        << "current 4vector:" << std::endl;
0069   for (int i = 0; i < 4; i++) {
0070     info << "p[" << i << "] = " << std::setw(20) << (*getIni4Vec())[i] << std::setw(20) << (*getCurr4Vec())[i]
0071          << std::endl;
0072   }
0073   info << "mass = " << std::setw(20) << (*getIni4Vec()).M() << std::setw(20) << (*getCurr4Vec()).M() << std::endl;
0074
0075   info << "u1  = " << _u1.X() << ", " << _u1.Y() << ", " << _u1.Z() << std::endl;
0076   info << "u2  = " << _u2.X() << ", " << _u2.Y() << ", " << _u2.Z() << std::endl;
0077   info << "u3  = " << _u3.X() << ", " << _u3.Y() << ", " << _u3.Z() << std::endl;
0078
0079   return info.str();
0080 }
0081
0082 void TAbsFitParticle::print() {
0083   // Print particle contents
0084
0085   edm::LogVerbatim("KinFitter") << this->getInfoString();
0086 }
0087
0088 void TAbsFitParticle::reset() {
0089   // Reset particle to initial values
0090
0091   _parameters = _iniparameters;
0092   _pcurr = _pini;
0093   setCovMatrixFit(nullptr);
0094   _pull.ResizeTo(_nPar, 1);
0095   _pull.Zero();
0096 }
0097
0098 void TAbsFitParticle::setCovMatrix(const TMatrixD* theCovMatrix) {
0099   // Set the measured covariance matrix
0100
0101   _covMatrix.ResizeTo(_nPar, _nPar);
0102   if (theCovMatrix == nullptr) {
0103     _covMatrix.Zero();
0104   } else if (theCovMatrix->GetNcols() == _nPar && theCovMatrix->GetNrows() == _nPar) {
0105     _covMatrix = (*theCovMatrix);
0106   } else {
0107     edm::LogError("WrongMatrixSize") << GetName() << "::setCovMatrix - Covariance matrix needs to be a " << _nPar << "x"
0108                                      << _nPar << " matrix.";
0109   }
0110 }
0111
0112 void TAbsFitParticle::setCovMatrixFit(const TMatrixD* theCovMatrixFit) {
0113   // Set the fitted covariance matrix
0114
0115   _covMatrixFit.ResizeTo(_nPar, _nPar);
0116   if (theCovMatrixFit == nullptr) {
0117     _covMatrixFit.Zero();
0118   } else if (theCovMatrixFit->GetNcols() == _nPar && theCovMatrixFit->GetNrows() == _nPar) {
0119     _covMatrixFit = (*theCovMatrixFit);
0120   } else {
0121     edm::LogError("WrongMatrixSize") << GetName() << "::setCovMatrixFit - Fitted covariance matrix needs to be a "
0122                                      << _nPar << "x" << _nPar << " matrix.";
0123   }
0124 }
0125
0126 void TAbsFitParticle::calcCovMatrixDeltaY() {
0127   // Calculates V(deltaY) ==  V(y_meas) - V(y_fit)
0128
0129   _covMatrixDeltaY.ResizeTo(_nPar, _nPar);
0130   _covMatrixDeltaY = _covMatrix;
0131   if (_covMatrixFit.GetNrows() == _nPar && _covMatrixFit.GetNcols() == _nPar)
0132     _covMatrixDeltaY -= _covMatrixFit;
0133   else
0134     edm::LogError("WrongMatrixSize") << GetName() << "::calcCovMatrixDeltaY - _covMatrixFit probably not set.";
0135 }
0136
0137 const TMatrixD* TAbsFitParticle::getPull() {
0138   // Calculates the pull (y_fit - y_meas) / sigma
0139   // with sigma = Sqrt( sigma[y_meas]^2 - V[y_fit]^2 )
0140   // for all parameters
0141
0142   _pull.ResizeTo(_nPar, 1);
0143   _pull = _parameters;
0144   _pull -= _iniparameters;
0145   calcCovMatrixDeltaY();
0146   for (int i = 0; i < _nPar; i++) {
0147     Double_t sigmaDeltaY = _covMatrixDeltaY(i, i);
0148     if (sigmaDeltaY < 0) {
0149       edm::LogWarning("NegativeDiagonalElem") << "V[deltaY] has a negative diagonal element.";
0150       _pull.Zero();
0151       return &_pull;
0152     } else {
0153       _pull(i, 0) /= TMath::Sqrt(sigmaDeltaY);
0154     }
0155   }
0156
0157   return &_pull;
0158 }
0159
0160 void TAbsFitParticle::applycorr(TMatrixD* corrMatrix) {
0161   // Apply corrections to the parameters wrt. to the
0162   // initial parameters y* = y + delta(y)
0163   // This method will also calculate the fitted
0164   // 4vector of the particle
0165
0166   // update _parameters-Matrix
0167   _parameters = _iniparameters;
0168   _parameters += (*corrMatrix);
0169
0170   // calculates new 4vec
0171   TLorentzVector* vec = calc4Vec(&_parameters);
0172   _pcurr = (*vec);
0173   delete vec;
0174 }
0175
0176 void TAbsFitParticle::setParIni(const TMatrixD* parini) {
0177   if (parini == nullptr)
0178     return;
0179   else if (parini->GetNrows() == _iniparameters.GetNrows() && parini->GetNcols() == _iniparameters.GetNcols())
0180     _iniparameters = (*parini);
0181   else {
0182     edm::LogError("WrongMatrixSize") << GetName() << "::setParIni - Matrices don't fit.";
0183     return;
0184   }
0185 }
0186
0187 const TMatrixD* TAbsFitParticle::getCovMatrixDeltaY() {
0188   //
0189   calcCovMatrixDeltaY();
0190   return &_covMatrixDeltaY;
0191 }