KinFitter/src/TFitParticleEScaledMomDev.cc

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