TAbsFitParticle.cc

Go to the documentation of this file.

// Classname: TAbsFitParticle
// Author: Jan E. Sundermann, Verena Klose (TU Dresden)      


//________________________________________________________________
// 
// TAbsFitParticle::
// --------------------
//
// Abstract base class for particles to be used with kinematic fitter
//

#include "PhysicsTools/KinFitter/interface/TAbsFitParticle.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include <iostream>
#include <iomanip>
#include "TClass.h"


TAbsFitParticle::TAbsFitParticle():
  TNamed("NoName","NoTitle")
  ,_nPar(0)     
  ,_u1()   
  ,_u2()
  ,_u3()
  ,_covMatrix()
  ,_covMatrixFit()
  ,_covMatrixDeltaY()
  ,_pull()
  ,_iniparameters(1,1)
  ,_parameters(1,1)
  ,_pini()
  ,_pcurr()
{
 
}

TAbsFitParticle::TAbsFitParticle(const TString &name, const TString &title ):
  TNamed(name,title)
  ,_nPar(0)     
  ,_u1()   
  ,_u2()
  ,_u3()
  ,_covMatrix()
  ,_covMatrixFit()
  ,_covMatrixDeltaY()
  ,_pull()
  ,_iniparameters(1,1)
  ,_parameters(1,1)
  ,_pini()
  ,_pcurr()
{
 
}


TAbsFitParticle::~TAbsFitParticle() {

}

TString
TAbsFitParticle::getInfoString() {
  // Collect information to be used for printout

  std::stringstream info;
  info << std::scientific << std::setprecision(6);

  info << "__________________________" << std::endl
       << std::endl;

  info << "OBJ: " << IsA()->GetName() << "\t" << GetName() << "\t" << GetTitle() << std::endl;

  info << std::setw(22) << "initial parameters:"  << std::setw(5) << " " << std::setw(20) << "current parameters:" << std::endl;
  for (int i = 0; i< _nPar ;i++){
    info << "par[" << i << "] = "
     << std::setw(18) << (*getParIni())(i,0) 
     << std::setw(20) << (*getParCurr())(i,0) << std::endl;
  }

  info << std::setw(22) << "initial 4vector:" << std::setw(5) << " " << std::setw(20) << "current 4vector:" << std::endl;
  for (int i = 0; i< 4 ;i++){
    info << "p[" << i << "] = "
     << std::setw(20) << (*getIni4Vec())[i] 
     << std::setw(20) << (*getCurr4Vec())[i] << std::endl;
   }
  info << "mass = " 
       << std::setw(20) << (*getIni4Vec()).M() 
       << std::setw(20) << (*getCurr4Vec()).M() << std::endl;

   info << "u1  = " << _u1.X() << ", " << _u1.Y() << ", " << _u1.Z() << std::endl;
   info << "u2  = " << _u2.X() << ", " << _u2.Y() << ", " << _u2.Z() << std::endl;
   info << "u3  = " << _u3.X() << ", " << _u3.Y() << ", " << _u3.Z() << std::endl;

   return info.str();

}

void 
TAbsFitParticle::print() {
  // Print particle contents

  edm::LogVerbatim("KinFitter") << this->getInfoString();

}

void TAbsFitParticle::reset() {
  // Reset particle to initial values

  _parameters = _iniparameters;  
  _pcurr = _pini;
  setCovMatrixFit( nullptr );
  _pull.ResizeTo(_nPar, 1);
  _pull.Zero();

}

void TAbsFitParticle::setCovMatrix(const TMatrixD* theCovMatrix) {
  // Set the measured covariance matrix

  _covMatrix.ResizeTo(_nPar, _nPar);
  if(theCovMatrix==nullptr) {
    _covMatrix.Zero();
  } else if (theCovMatrix->GetNcols() ==_nPar && theCovMatrix->GetNrows() ==_nPar) {
    _covMatrix = (*theCovMatrix);
  } else {
    edm::LogError ("WrongMatrixSize")
      << GetName() << "::setCovMatrix - Covariance matrix needs to be a "
      << _nPar << "x" << _nPar << " matrix.";
  }

}


void TAbsFitParticle::setCovMatrixFit(const TMatrixD* theCovMatrixFit) {
  // Set the fitted covariance matrix

  _covMatrixFit.ResizeTo(_nPar, _nPar);
  if(theCovMatrixFit==nullptr) {
    _covMatrixFit.Zero();
  } else if (theCovMatrixFit->GetNcols() ==_nPar && theCovMatrixFit->GetNrows() ==_nPar) {
    _covMatrixFit = (*theCovMatrixFit);
  } else {
    edm::LogError ("WrongMatrixSize")
      << GetName() << "::setCovMatrixFit - Fitted covariance matrix needs to be a "
      << _nPar << "x" << _nPar << " matrix.";
  }

}

void TAbsFitParticle::calcCovMatrixDeltaY() {
  // Calculates V(deltaY) ==  V(y_meas) - V(y_fit)

  _covMatrixDeltaY.ResizeTo( _nPar, _nPar );
  _covMatrixDeltaY = _covMatrix;
  if(_covMatrixFit.GetNrows() == _nPar && _covMatrixFit.GetNcols() == _nPar)
    _covMatrixDeltaY -= _covMatrixFit;
  else 
    edm::LogError ("WrongMatrixSize")
      << GetName() << "::calcCovMatrixDeltaY - _covMatrixFit probably not set.";
}

const TMatrixD* TAbsFitParticle::getPull() {
  // Calculates the pull (y_fit - y_meas) / sigma
  // with sigma = Sqrt( sigma[y_meas]^2 - V[y_fit]^2 )
  // for all parameters

  _pull.ResizeTo( _nPar, 1 );
  _pull = _parameters;
  _pull -= _iniparameters;
  calcCovMatrixDeltaY(); 
  for (int i = 0; i<_nPar; i++) {
    Double_t sigmaDeltaY = _covMatrixDeltaY(i, i);
    if (sigmaDeltaY < 0) {
      edm::LogWarning ("NegativeDiagonalElem") << "V[deltaY] has a negative diagonal element.";
      _pull.Zero();
      return &_pull;
    } else {
      _pull(i,0) /= TMath::Sqrt( sigmaDeltaY );
    }
  }

  return &_pull;

}

void TAbsFitParticle::applycorr(TMatrixD* corrMatrix) {
  // Apply corrections to the parameters wrt. to the
  // initial parameters y* = y + delta(y)
  // This method will also calculate the fitted 
  // 4vector of the particle

  // update _parameters-Matrix
  _parameters = _iniparameters;
  _parameters += (*corrMatrix);

  // calculates new 4vec
  TLorentzVector* vec = calc4Vec( &_parameters );
  _pcurr = (*vec);
  delete vec;

}

void TAbsFitParticle::setParIni(const TMatrixD* parini) {
  if (parini == nullptr) return;
  else if( parini->GetNrows() == _iniparameters.GetNrows()
       && parini->GetNcols() == _iniparameters.GetNcols() )
    _iniparameters = (*parini) ;
  else {
    edm::LogError ("WrongMatrixSize")
      << GetName() << "::setParIni - Matrices don't fit.";
    return;
      }
}

const TMatrixD* TAbsFitParticle::getCovMatrixDeltaY() {
  //
  calcCovMatrixDeltaY(); 
  return &_covMatrixDeltaY; 
}