KinFitter/src/TFitParticleEMomDev.cc

0001 // Classname: TFitParticleEMomDev
0002 // Author: Jan E. Sundermann, Verena Klose (TU Dresden)
0003
0004 //________________________________________________________________
0005 //
0006 // TFitParticleEMomDev
0007 // --------------------
0008 //
0009 // Particle with special parametrization of the momentum 4vector and
0010 // free mass (4 free parameters). The parametrization is chosen as
0011 // follows:
0012 //
0013 // p = r*|p|*u_r + theta*u_theta + phi*u_phi
0014 // E(fit) = E_meas * d
0015 //
0016 // with u_r = p/|p|
0017 //      u_phi = (u_z x u_r)/|u_z x u_r|
0018 //      u_theta = (u_r x u_phi)/|u_r x u_phi|
0019 //
0020 // The initial parameters values are chosen like (r, theta, phi, d) = (1., 0., 0., 1.)
0021 // corresponding to the measured momentum and mass.
0022 //
0023
0024 #include <iostream>
0025 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0026 #include "PhysicsTools/KinFitter/interface/TFitParticleEMomDev.h"
0027
0028 //----------------
0029 // Constructor --
0030 //----------------
0031 TFitParticleEMomDev::TFitParticleEMomDev() : TAbsFitParticle() { init(nullptr, nullptr); }
0032
0033 TFitParticleEMomDev::TFitParticleEMomDev(const TFitParticleEMomDev& fitParticle)
0034     : TAbsFitParticle(fitParticle.GetName(), fitParticle.GetTitle()) {
0035   _nPar = fitParticle._nPar;
0036   _u1 = fitParticle._u1;
0037   _u2 = fitParticle._u2;
0038   _u3 = fitParticle._u3;
0039   _covMatrix.ResizeTo(fitParticle._covMatrix);
0040   _covMatrix = fitParticle._covMatrix;
0041   _iniparameters.ResizeTo(fitParticle._iniparameters);
0042   _iniparameters = fitParticle._iniparameters;
0043   _parameters.ResizeTo(fitParticle._parameters);
0044   _parameters = fitParticle._parameters;
0045   _pini = fitParticle._pini;
0046   _pcurr = fitParticle._pcurr;
0047 }
0048
0049 TFitParticleEMomDev::TFitParticleEMomDev(TLorentzVector* pini, const TMatrixD* theCovMatrix) : TAbsFitParticle() {
0050   init(pini, theCovMatrix);
0051 }
0052
0053 TFitParticleEMomDev::TFitParticleEMomDev(const TString& name,
0054                                          const TString& title,
0055                                          TLorentzVector* pini,
0056                                          const TMatrixD* theCovMatrix)
0057     : TAbsFitParticle(name, title) {
0058   init(pini, theCovMatrix);
0059 }
0060
0061 TAbsFitParticle* TFitParticleEMomDev::clone(const TString& newname) const {
0062   // Returns a copy of itself
0063
0064   TAbsFitParticle* myclone = new TFitParticleEMomDev(*this);
0065   if (newname.Length() > 0)
0066     myclone->SetName(newname);
0067   return myclone;
0068 }
0069
0070 //--------------
0071 // Destructor --
0072 //--------------
0073 TFitParticleEMomDev::~TFitParticleEMomDev() {}
0074
0075 //--------------
0076 // Operations --
0077 //--------------
0078 void TFitParticleEMomDev::init(TLorentzVector* pini, const TMatrixD* theCovMatrix) {
0079   _nPar = 4;
0080   setIni4Vec(pini);
0081   setCovMatrix(theCovMatrix);
0082 }
0083
0084 void TFitParticleEMomDev::setIni4Vec(const TLorentzVector* pini) {
0085   // Set the initial 4vector. Will also set the
0086   // inital parameter values
0087   _iniparameters.ResizeTo(_nPar, 1);
0088   _parameters.ResizeTo(_nPar, 1);
0089
0090   if (pini == nullptr) {
0091     _iniparameters(0, 0) = 0.;
0092     _iniparameters(1, 0) = 0.;
0093     _iniparameters(2, 0) = 0.;
0094     _iniparameters(3, 0) = 0.;
0095
0096     _u1.SetXYZ(0., 0., 0.);
0097     _u3.SetXYZ(0., 0., 0.);
0098     _u2.SetXYZ(0., 0., 0.);
0099     _pini.SetXYZT(0., 0., 0., 0.);
0100     _pcurr = _pini;
0101
0102   } else {
0103     _iniparameters(0, 0) = 1.;
0104     _iniparameters(1, 0) = 0.;
0105     _iniparameters(2, 0) = 0.;
0106     _iniparameters(3, 0) = 1.;
0107
0108     _pini = (*pini);
0109     _pcurr = _pini;
0110
0111     _u1 = pini->Vect();
0112     _u1 *= 1. / _u1.Mag();
0113
0114     TVector3 uz(0., 0., 1.);
0115     _u3 = uz.Cross(_u1);
0116     _u3 *= 1. / _u3.Mag();
0117
0118     _u2 = _u3.Cross(_u1);
0119     _u2 *= 1. / _u2.Mag();
0120   }
0121
0122   // reset parameters
0123   _parameters = _iniparameters;
0124 }
0125
0126 TLorentzVector* TFitParticleEMomDev::calc4Vec(const TMatrixD* params) {
0127   // Calculates a 4vector corresponding to the given
0128   // parameter values
0129
0130   if (params == nullptr) {
0131     return nullptr;
0132   }
0133
0134   if (params->GetNcols() != 1 || params->GetNrows() != _nPar) {
0135     edm::LogError("WrongMatrixSize") << GetName() << "::calc4Vec - Parameter matrix has wrong size.";
0136     return nullptr;
0137   }
0138
0139   Double_t X = _pini.P() * (*params)(0, 0) * _u1.X() + (*params)(1, 0) * _u2.X() + (*params)(2, 0) * _u3.X();
0140   Double_t Y = _pini.P() * (*params)(0, 0) * _u1.Y() + (*params)(1, 0) * _u2.Y() + (*params)(2, 0) * _u3.Y();
0141   Double_t Z = _pini.P() * (*params)(0, 0) * _u1.Z() + (*params)(1, 0) * _u2.Z() + (*params)(2, 0) * _u3.Z();
0142   Double_t E = _pini.E() * (*params)(3, 0);
0143
0144   TLorentzVector* vec = new TLorentzVector(X, Y, Z, E);
0145   return vec;
0146 }
0147
0148 TMatrixD* TFitParticleEMomDev::getDerivative() {
0149   // returns derivative dP/dy with P=(p,E) and y=(r, theta, phi, d)
0150   // the free parameters of the fit. The columns of the matrix contain
0151   // (dP/dr, dP/dtheta, dP/dphi, dP/dd).
0152
0153   TMatrixD* DerivativeMatrix = new TMatrixD(4, 4);
0154   (*DerivativeMatrix) *= 0.;
0155
0156   //1st column: dP/dr
0157   (*DerivativeMatrix)(0, 0) = _pini.P() * _u1.X();
0158   (*DerivativeMatrix)(1, 0) = _pini.P() * _u1.Y();
0159   (*DerivativeMatrix)(2, 0) = _pini.P() * _u1.Z();
0160   (*DerivativeMatrix)(3, 0) = 0.;
0161
0162   //2nd column: dP/dtheta
0163   (*DerivativeMatrix)(0, 1) = _u2.X();
0164   (*DerivativeMatrix)(1, 1) = _u2.Y();
0165   (*DerivativeMatrix)(2, 1) = _u2.Z();
0166   (*DerivativeMatrix)(3, 1) = 0.;
0167
0168   //3rd column: dP/dphi
0169   (*DerivativeMatrix)(0, 2) = _u3.X();
0170   (*DerivativeMatrix)(1, 2) = _u3.Y();
0171   (*DerivativeMatrix)(2, 2) = _u3.Z();
0172   (*DerivativeMatrix)(3, 2) = 0.;
0173
0174   //4th column: dP/dm
0175   (*DerivativeMatrix)(0, 3) = 0.;
0176   (*DerivativeMatrix)(1, 3) = 0.;
0177   (*DerivativeMatrix)(2, 3) = 0.;
0178   (*DerivativeMatrix)(3, 3) = _pini.E();
0179
0180   return DerivativeMatrix;
0181 }
0182
0183 TMatrixD* TFitParticleEMomDev::transform(const TLorentzVector& vec) {
0184   // Returns the parameters corresponding to the given
0185   // 4vector wrt. to the current base vectors u_r, u_theta, and u_phi
0186
0187   // construct rotation matrix
0188   TRotation rot;
0189   rot.RotateAxes(_u1, _u2, _u3);
0190   rot.Invert();
0191
0192   // rotate vector
0193   TVector3 vec3(vec.Vect());
0194   vec3.Transform(rot);
0195
0196   // retrieve parameters
0197   TMatrixD* tparams = new TMatrixD(_nPar, 1);
0198   (*tparams)(0, 0) = vec3.x() / _pini.P();
0199   (*tparams)(1, 0) = vec3.y();
0200   (*tparams)(2, 0) = vec3.z();
0201   (*tparams)(3, 0) = vec.E() / _pini.E();
0202
0203   return tparams;
0204 }