KinFitter/src/TFitParticleMCSpher.cc

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