KinFitter/src/TFitParticleMomDev.cc

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