KinFitter/src/TFitParticleMCMomDev.cc

0001 // Classname: TFitParticleMCMomDev
0002 // Author: Jan E. Sundermann, Verena Klose (TU Dresden)
0003
0004 //________________________________________________________________
0005 //
0006 // TFitParticleMCMomDev
0007 // --------------------
0008 //
0009 // Particle with special parametrization of the momentum 4vector and
0010 // constant mass (3 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 + 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) = (1., 0., 0.)
0021 // corresponding to the measured momentum.
0022 //
0023
0024 #include <iostream>
0025 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0026 #include "PhysicsTools/KinFitter/interface/TFitParticleMCMomDev.h"
0027 #include "TMath.h"
0028
0029 //----------------
0030 // Constructor --
0031 //----------------
0032 TFitParticleMCMomDev::TFitParticleMCMomDev() : TAbsFitParticle() { init(nullptr, 0., nullptr); }
0033
0034 TFitParticleMCMomDev::TFitParticleMCMomDev(const TFitParticleMCMomDev& 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 TFitParticleMCMomDev::TFitParticleMCMomDev(TVector3* p, Double_t M, const TMatrixD* theCovMatrix) : TAbsFitParticle() {
0051   init(p, M, theCovMatrix);
0052 }
0053
0054 TFitParticleMCMomDev::TFitParticleMCMomDev(
0055     const TString& name, const TString& title, TVector3* p, Double_t M, const TMatrixD* theCovMatrix)
0056     : TAbsFitParticle(name, title) {
0057   init(p, M, theCovMatrix);
0058 }
0059
0060 TAbsFitParticle* TFitParticleMCMomDev::clone(const TString& newname) const {
0061   // Returns a copy of itself
0062
0063   TAbsFitParticle* myclone = new TFitParticleMCMomDev(*this);
0064   if (newname.Length() > 0)
0065     myclone->SetName(newname);
0066   return myclone;
0067 }
0068
0069 //--------------
0070 // Destructor --
0071 //--------------
0072 TFitParticleMCMomDev::~TFitParticleMCMomDev() {}
0073
0074 //--------------
0075 // Operations --
0076 //--------------
0077 void TFitParticleMCMomDev::init(TVector3* p, Double_t M, const TMatrixD* theCovMatrix) {
0078   _nPar = 3;
0079   setIni4Vec(p, M);
0080   _iniparameters.ResizeTo(_nPar, 1);
0081   _iniparameters(0, 0) = 1.;
0082   _iniparameters(1, 0) = 0.;
0083   _iniparameters(2, 0) = 0.;
0084   _parameters.ResizeTo(_nPar, 1);
0085   _parameters = _iniparameters;
0086   setCovMatrix(theCovMatrix);
0087 }
0088
0089 TLorentzVector* TFitParticleMCMomDev::calc4Vec(const TMatrixD* params) {
0090   // Calculates a 4vector corresponding to the given
0091   // parameter values
0092
0093   if (params == nullptr) {
0094     return nullptr;
0095   }
0096
0097   if (params->GetNcols() != 1 || params->GetNrows() != _nPar) {
0098     edm::LogError("WrongMatrixSize") << GetName() << "::calc4Vec - Parameter matrix has wrong size.";
0099     return nullptr;
0100   }
0101
0102   Double_t X = _pini.P() * (*params)(0, 0) * _u1.X() + (*params)(1, 0) * _u2.X() + (*params)(2, 0) * _u3.X();
0103   Double_t Y = _pini.P() * (*params)(0, 0) * _u1.Y() + (*params)(1, 0) * _u2.Y() + (*params)(2, 0) * _u3.Y();
0104   Double_t Z = _pini.P() * (*params)(0, 0) * _u1.Z() + (*params)(1, 0) * _u2.Z() + (*params)(2, 0) * _u3.Z();
0105   Double_t E = TMath::Sqrt(X * X + Y * Y + Z * Z + _pini.M2());
0106
0107   TLorentzVector* vec = new TLorentzVector(X, Y, Z, E);
0108   return vec;
0109 }
0110
0111 void TFitParticleMCMomDev::setIni4Vec(const TLorentzVector* pini) {
0112   // Set the initial 4vector. Will also set the
0113   // inital parameter values
0114
0115   TVector3 vec(pini->Vect());
0116   setIni4Vec(&vec, pini->M());
0117 }
0118
0119 void TFitParticleMCMomDev::setIni4Vec(const TVector3* p, Double_t M) {
0120   // Set the initial 4vector. Will also set the
0121   // inital parameter values
0122
0123   if (p == nullptr) {
0124     _pini.SetXYZM(0., 0., 0., M);
0125     _pcurr = _pini;
0126
0127   } else {
0128     _pini.SetXYZM(p->x(), p->y(), p->z(), M);
0129     _pcurr = _pini;
0130
0131     _u1 = (*p);
0132     _u1 *= 1. / _u1.Mag();
0133
0134     TVector3 uz(0., 0., 1.);
0135     _u3 = uz.Cross(_u1);
0136     _u3 *= 1. / _u3.Mag();
0137
0138     _u2 = _u3.Cross(_u1);
0139     _u2 *= 1. / _u2.Mag();
0140   }
0141
0142   _parameters = _iniparameters;
0143 }
0144
0145 TMatrixD* TFitParticleMCMomDev::getDerivative() {
0146   // returns derivative dP/dy with P=(p,E) and y=(r, theta, phi)
0147   // the free parameters of the fit. The columns of the matrix contain
0148   // (dP/dr, dP/dtheta, dP/dphi).
0149
0150   TMatrixD* DerivativeMatrix = new TMatrixD(4, 3);
0151   (*DerivativeMatrix) *= 0.;
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) = _pini.P() * _pini.P() * _parameters(0, 0) / _pcurr.E();
0157
0158   //  (*DerivativeMatrix)(3,0)=0.;
0159
0160   //2nd column: dP/dtheta
0161   (*DerivativeMatrix)(0, 1) = _u2.X();
0162   (*DerivativeMatrix)(1, 1) = _u2.Y();
0163   (*DerivativeMatrix)(2, 1) = _u2.Z();
0164   (*DerivativeMatrix)(3, 1) = _parameters(1, 0) / _pcurr.E();
0165
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) = _parameters(2, 0) / _pcurr.E();
0173
0174   //(*DerivativeMatrix)(3,2)=0.;
0175
0176   return DerivativeMatrix;
0177 }
0178
0179 TMatrixD* TFitParticleMCMomDev::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
0198   return tparams;
0199 }