KinFitter/src/TFitParticleMCPInvSpher.cc

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