core/src/AbsTrackRep.cc

0001 /* Copyright 2008-2009, Technische Universitaet Muenchen,
0002    Authors: Christian Hoeppner & Sebastian Neubert
0003
0004    This file is part of GENFIT.
0005
0006    GENFIT is free software: you can redistribute it and/or modify
0007    it under the terms of the GNU Lesser General Public License as published
0008    by the Free Software Foundation, either version 3 of the License, or
0009    (at your option) any later version.
0010
0011    GENFIT is distributed in the hope that it will be useful,
0012    but WITHOUT ANY WARRANTY; without even the implied warranty of
0013    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0014    GNU Lesser General Public License for more details.
0015
0016    You should have received a copy of the GNU Lesser General Public License
0017    along with GENFIT.  If not, see <http://www.gnu.org/licenses/>.
0018 */
0019
0020 #include "AbsTrackRep.h"
0021 #include "StateOnPlane.h"
0022 #include "AbsMeasurement.h"
0023 #include "IO.h"
0024
0025 #include <TDatabasePDG.h>
0026
0027
0028
0029 namespace genfit {
0030
0031 AbsTrackRep::AbsTrackRep() :
0032   pdgCode_(0), propDir_(0), debugLvl_(0)
0033 {
0034   ;
0035 }
0036
0037 AbsTrackRep::AbsTrackRep(int pdgCode, char propDir) :
0038   pdgCode_(pdgCode), propDir_(propDir), debugLvl_(0)
0039 {
0040   ;
0041 }
0042
0043 AbsTrackRep::AbsTrackRep(const AbsTrackRep& rep) :
0044   TObject(rep), pdgCode_(rep.pdgCode_), propDir_(rep.propDir_), debugLvl_(rep.debugLvl_)
0045 {
0046   ;
0047 }
0048
0049
0050 double AbsTrackRep::extrapolateToMeasurement(StateOnPlane& state,
0051     const AbsMeasurement* measurement,
0052     bool stopAtBoundary,
0053     bool calcJacobianNoise) const {
0054
0055   return this->extrapolateToPlane(state, measurement->constructPlane(state), stopAtBoundary, calcJacobianNoise);
0056 }
0057
0058
0059 TVectorD AbsTrackRep::get6DState(const StateOnPlane& state) const {
0060   TVector3 pos, mom;
0061   getPosMom(state, pos, mom);
0062
0063   TVectorD stateVec(6);
0064
0065   stateVec(0) = pos.X();
0066   stateVec(1) = pos.Y();
0067   stateVec(2) = pos.Z();
0068
0069   stateVec(3) = mom.X();
0070   stateVec(4) = mom.Y();
0071   stateVec(5) = mom.Z();
0072
0073   return stateVec;
0074 }
0075
0076
0077 void AbsTrackRep::get6DStateCov(const MeasuredStateOnPlane& state, TVectorD& stateVec, TMatrixDSym& cov) const {
0078   TVector3 pos, mom;
0079   getPosMomCov(state, pos, mom, cov);
0080
0081   stateVec.ResizeTo(6);
0082
0083   stateVec(0) = pos.X();
0084   stateVec(1) = pos.Y();
0085   stateVec(2) = pos.Z();
0086
0087   stateVec(3) = mom.X();
0088   stateVec(4) = mom.Y();
0089   stateVec(5) = mom.Z();
0090 }
0091
0092
0093 double AbsTrackRep::getPDGCharge() const {
0094   TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(pdgCode_);
0095   assert(particle != nullptr);
0096   return particle->Charge()/(3.);
0097 }
0098
0099
0100 double AbsTrackRep::getMass(const StateOnPlane& /*state*/) const {
0101   return TDatabasePDG::Instance()->GetParticle(pdgCode_)->Mass();
0102 }
0103
0104
0105 void AbsTrackRep::calcJacobianNumerically(const genfit::StateOnPlane& origState,
0106                                                const genfit::SharedPlanePtr destPlane,
0107                                                TMatrixD& jacobian) const {
0108
0109   // Find the transport matrix for track propagation from origState to destPlane
0110   // I.e. this finds
0111   //            d stateDestPlane / d origState |_origState
0112
0113   jacobian.ResizeTo(getDim(), getDim());
0114
0115   // no science behind these values, I verified that forward and
0116   // backward propagation yield inverse matrices to good
0117   // approximation.  In order to avoid bad roundoff errors, the actual
0118   // step taken is determined below, separately for each direction.
0119   const double defaultStepX = 1.E-5;
0120   double stepX;
0121
0122   // Calculate derivative for all three dimensions successively.
0123   // The algorithm follows the one in TF1::Derivative() :
0124   //   df(x) = (4 D(h/2) - D(h)) / 3
0125   // with D(h) = (f(x + h) - f(x - h)) / (2 h).
0126   //
0127   // Could perhaps do better by also using f(x) which would be stB.
0128   TVectorD rightShort(getDim()), rightFull(getDim());
0129   TVectorD leftShort(getDim()), leftFull(getDim());
0130   for (size_t i = 0; i < getDim(); ++i) {
0131     {
0132       genfit::StateOnPlane stateCopy(origState);
0133       double temp = stateCopy.getState()(i) + defaultStepX / 2;
0134       // Find the actual size of the step, which will differ from
0135       // defaultStepX due to roundoff.  This is the step-size we will
0136       // use for this direction.  Idea taken from Numerical Recipes,
0137       // 3rd ed., section 5.7.
0138       //
0139       // Note that if a number is exactly representable, it's double
0140       // will also be exact.  Outside denormals, this also holds for
0141       // halving.  Unless the exponent changes (which it only will in
0142       // the vicinity of zero) adding or subtracing doesn't make a
0143       // difference.
0144       //
0145       // We determine the roundoff error for the half-step.  If this
0146       // is exactly representable, the full step will also be.
0147       stepX = 2 * (temp - stateCopy.getState()(i));
0148       (stateCopy.getState())(i) = temp;
0149       extrapolateToPlane(stateCopy, destPlane);
0150       rightShort = stateCopy.getState();
0151     }
0152     {
0153       genfit::StateOnPlane stateCopy(origState);
0154       (stateCopy.getState())(i) -= stepX / 2;
0155       extrapolateToPlane(stateCopy, destPlane);
0156       leftShort = stateCopy.getState();
0157     }
0158     {
0159       genfit::StateOnPlane stateCopy(origState);
0160       (stateCopy.getState())(i) += stepX;
0161       extrapolateToPlane(stateCopy, destPlane);
0162       rightFull = stateCopy.getState();
0163     }
0164     {
0165       genfit::StateOnPlane stateCopy(origState);
0166       (stateCopy.getState())(i) -= stepX;
0167       extrapolateToPlane(stateCopy, destPlane);
0168       leftFull = stateCopy.getState();
0169     }
0170
0171     // Calculate the derivatives for the individual components of
0172     // the track parameters.
0173     for (size_t j = 0; j < getDim(); ++j) {
0174       double derivFull = (rightFull(j) - leftFull(j)) / 2 / stepX;
0175       double derivShort = (rightShort(j) - leftShort(j)) / stepX;
0176
0177       jacobian(j, i) = 1./3.*(4*derivShort - derivFull);
0178     }
0179   }
0180 }
0181
0182
0183 bool AbsTrackRep::switchPDGSign() {
0184   TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(-pdgCode_);
0185   if(particle != nullptr) {
0186     pdgCode_ *= -1;
0187     return true;
0188   }
0189   return false;
0190 }
0191
0192
0193
0194 void AbsTrackRep::Print(const Option_t*) const {
0195   printOut << "genfit::AbsTrackRep, pdgCode = " << pdgCode_ << ". PropDir = " << (int)propDir_ << "\n";
0196 }
0197
0198
0199 } /* End of namespace genfit */