Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AbsTrackRep.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file AbsTrackRep.cc
1 /* Copyright 2008-2009, Technische Universitaet Muenchen,
2  Authors: Christian Hoeppner & Sebastian Neubert
3 
4  This file is part of GENFIT.
5 
6  GENFIT is free software: you can redistribute it and/or modify
7  it under the terms of the GNU Lesser General Public License as published
8  by the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  GENFIT is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU Lesser General Public License for more details.
15 
16  You should have received a copy of the GNU Lesser General Public License
17  along with GENFIT. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 #include "AbsTrackRep.h"
21 #include "StateOnPlane.h"
22 #include "AbsMeasurement.h"
23 #include "IO.h"
24 
25 #include <TDatabasePDG.h>
26 
27 
28 
29 namespace genfit {
30 
32  pdgCode_(0), propDir_(0), debugLvl_(0)
33 {
34  ;
35 }
36 
37 AbsTrackRep::AbsTrackRep(int pdgCode, char propDir) :
38  pdgCode_(pdgCode), propDir_(propDir), debugLvl_(0)
39 {
40  ;
41 }
42 
44  TObject(rep), pdgCode_(rep.pdgCode_), propDir_(rep.propDir_), debugLvl_(rep.debugLvl_)
45 {
46  ;
47 }
48 
49 
51  const AbsMeasurement* measurement,
52  bool stopAtBoundary,
53  bool calcJacobianNoise) const {
54 
55  return this->extrapolateToPlane(state, measurement->constructPlane(state), stopAtBoundary, calcJacobianNoise);
56 }
57 
58 
59 TVectorD AbsTrackRep::get6DState(const StateOnPlane& state) const {
60  TVector3 pos, mom;
61  getPosMom(state, pos, mom);
62 
63  TVectorD stateVec(6);
64 
65  stateVec(0) = pos.X();
66  stateVec(1) = pos.Y();
67  stateVec(2) = pos.Z();
68 
69  stateVec(3) = mom.X();
70  stateVec(4) = mom.Y();
71  stateVec(5) = mom.Z();
72 
73  return stateVec;
74 }
75 
76 
77 void AbsTrackRep::get6DStateCov(const MeasuredStateOnPlane& state, TVectorD& stateVec, TMatrixDSym& cov) const {
78  TVector3 pos, mom;
79  getPosMomCov(state, pos, mom, cov);
80 
81  stateVec.ResizeTo(6);
82 
83  stateVec(0) = pos.X();
84  stateVec(1) = pos.Y();
85  stateVec(2) = pos.Z();
86 
87  stateVec(3) = mom.X();
88  stateVec(4) = mom.Y();
89  stateVec(5) = mom.Z();
90 }
91 
92 
93 double AbsTrackRep::getPDGCharge() const {
94  TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(pdgCode_);
95  assert(particle != nullptr);
96  return particle->Charge()/(3.);
97 }
98 
99 
100 double AbsTrackRep::getMass(const StateOnPlane& /*state*/) const {
101  return TDatabasePDG::Instance()->GetParticle(pdgCode_)->Mass();
102 }
103 
104 
106  const genfit::SharedPlanePtr destPlane,
107  TMatrixD& jacobian) const {
108 
109  // Find the transport matrix for track propagation from origState to destPlane
110  // I.e. this finds
111  // d stateDestPlane / d origState |_origState
112 
113  jacobian.ResizeTo(getDim(), getDim());
114 
115  // no science behind these values, I verified that forward and
116  // backward propagation yield inverse matrices to good
117  // approximation. In order to avoid bad roundoff errors, the actual
118  // step taken is determined below, separately for each direction.
119  const double defaultStepX = 1.E-5;
120  double stepX;
121 
122  // Calculate derivative for all three dimensions successively.
123  // The algorithm follows the one in TF1::Derivative() :
124  // df(x) = (4 D(h/2) - D(h)) / 3
125  // with D(h) = (f(x + h) - f(x - h)) / (2 h).
126  //
127  // Could perhaps do better by also using f(x) which would be stB.
128  TVectorD rightShort(getDim()), rightFull(getDim());
129  TVectorD leftShort(getDim()), leftFull(getDim());
130  for (size_t i = 0; i < getDim(); ++i) {
131  {
132  genfit::StateOnPlane stateCopy(origState);
133  double temp = stateCopy.getState()(i) + defaultStepX / 2;
134  // Find the actual size of the step, which will differ from
135  // defaultStepX due to roundoff. This is the step-size we will
136  // use for this direction. Idea taken from Numerical Recipes,
137  // 3rd ed., section 5.7.
138  //
139  // Note that if a number is exactly representable, it's double
140  // will also be exact. Outside denormals, this also holds for
141  // halving. Unless the exponent changes (which it only will in
142  // the vicinity of zero) adding or subtracing doesn't make a
143  // difference.
144  //
145  // We determine the roundoff error for the half-step. If this
146  // is exactly representable, the full step will also be.
147  stepX = 2 * (temp - stateCopy.getState()(i));
148  (stateCopy.getState())(i) = temp;
149  extrapolateToPlane(stateCopy, destPlane);
150  rightShort = stateCopy.getState();
151  }
152  {
153  genfit::StateOnPlane stateCopy(origState);
154  (stateCopy.getState())(i) -= stepX / 2;
155  extrapolateToPlane(stateCopy, destPlane);
156  leftShort = stateCopy.getState();
157  }
158  {
159  genfit::StateOnPlane stateCopy(origState);
160  (stateCopy.getState())(i) += stepX;
161  extrapolateToPlane(stateCopy, destPlane);
162  rightFull = stateCopy.getState();
163  }
164  {
165  genfit::StateOnPlane stateCopy(origState);
166  (stateCopy.getState())(i) -= stepX;
167  extrapolateToPlane(stateCopy, destPlane);
168  leftFull = stateCopy.getState();
169  }
170 
171  // Calculate the derivatives for the individual components of
172  // the track parameters.
173  for (size_t j = 0; j < getDim(); ++j) {
174  double derivFull = (rightFull(j) - leftFull(j)) / 2 / stepX;
175  double derivShort = (rightShort(j) - leftShort(j)) / stepX;
176 
177  jacobian(j, i) = 1./3.*(4*derivShort - derivFull);
178  }
179  }
180 }
181 
182 
184  TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(-pdgCode_);
185  if(particle != nullptr) {
186  pdgCode_ *= -1;
187  return true;
188  }
189  return false;
190 }
191 
192 
193 
194 void AbsTrackRep::Print(const Option_t*) const {
195  printOut << "genfit::AbsTrackRep, pdgCode = " << pdgCode_ << ". PropDir = " << (int)propDir_ << "\n";
196 }
197 
198 
199 } /* End of namespace genfit */