Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RKTrackRep.h
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RKTrackRep.h
1 /* Copyright 2008-2010, Technische Universitaet Muenchen,
2  Authors: Christian Hoeppner & Sebastian Neubert & Johannes Rauch
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 
24 #ifndef genfit_RKTrackRep_h
25 #define genfit_RKTrackRep_h
26 
27 #include "AbsTrackRep.h"
28 #include "StateOnPlane.h"
29 #include "RKTools.h"
30 #include "StepLimits.h"
31 #include "Material.h"
32 
33 #include <algorithm>
34 
35 namespace genfit {
36 
40 struct RKStep {
41  MatStep matStep_; // material properties and stepsize
42  M1x7 state7_; // 7D state vector
43  StepLimits limits_;
44 
45  RKStep() {
46  std::fill(state7_.begin(), state7_.end(), 0);
47  }
48 };
49 
50 
54 struct ExtrapStep {
55  M7x7 jac7_; // 5D jacobian of transport
56  M7x7 noise7_; // 5D noise matrix
57 
58  ExtrapStep() {
59  std::fill(jac7_.begin(), jac7_.end(), 0);
60  std::fill(noise7_.begin(), noise7_.end(), 0);
61  }
62 };
63 
64 
72 class RKTrackRep : public AbsTrackRep {
73  friend class RKTrackRepTests_momMag_Test;
74  friend class RKTrackRepTests_calcForwardJacobianAndNoise_Test;
75  friend class RKTrackRepTests_getState7_Test;
76  friend class RKTrackRepTests_getState5_Test;
77 
78  public:
79 
80  RKTrackRep();
81  RKTrackRep(int pdgCode, char propDir = 0);
82 
83  virtual ~RKTrackRep();
84 
85  virtual AbsTrackRep* clone() const override {return new RKTrackRep(*this);}
86 
87  virtual double extrapolateToPlane(StateOnPlane& state,
88  const SharedPlanePtr& plane,
89  bool stopAtBoundary = false,
90  bool calcJacobianNoise = false) const override;
91 
92  using AbsTrackRep::extrapolateToLine;
93 
94  virtual double extrapolateToLine(StateOnPlane& state,
95  const TVector3& linePoint,
96  const TVector3& lineDirection,
97  bool stopAtBoundary = false,
98  bool calcJacobianNoise = false) const override;
99 
100  virtual double extrapolateToPoint(StateOnPlane& state,
101  const TVector3& point,
102  bool stopAtBoundary = false,
103  bool calcJacobianNoise = false) const override {
104  return extrapToPoint(state, point, nullptr, stopAtBoundary, calcJacobianNoise);
105  }
106 
107  virtual double extrapolateToPoint(StateOnPlane& state,
108  const TVector3& point,
109  const TMatrixDSym& G, // weight matrix (metric)
110  bool stopAtBoundary = false,
111  bool calcJacobianNoise = false) const override {
112  return extrapToPoint(state, point, &G, stopAtBoundary, calcJacobianNoise);
113  }
114 
115  virtual double extrapolateToCylinder(StateOnPlane& state,
116  double radius,
117  const TVector3& linePoint = TVector3(0.,0.,0.),
118  const TVector3& lineDirection = TVector3(0.,0.,1.),
119  bool stopAtBoundary = false,
120  bool calcJacobianNoise = false) const override;
121 
122 
123  virtual double extrapolateToCone(StateOnPlane& state,
124  double radius,
125  const TVector3& linePoint = TVector3(0.,0.,0.),
126  const TVector3& lineDirection = TVector3(0.,0.,1.),
127  bool stopAtBoundary = false,
128  bool calcJacobianNoise = false) const override ;
129 
130  virtual double extrapolateToSphere(StateOnPlane& state,
131  double radius,
132  const TVector3& point = TVector3(0.,0.,0.),
133  bool stopAtBoundary = false,
134  bool calcJacobianNoise = false) const override;
135 
136  virtual double extrapolateBy(StateOnPlane& state,
137  double step,
138  bool stopAtBoundary = false,
139  bool calcJacobianNoise = false) const override;
140 
141 
142  unsigned int getDim() const override {return 5;}
143 
144  virtual TVector3 getPos(const StateOnPlane& state) const override;
145 
146  virtual TVector3 getMom(const StateOnPlane& state) const override;
147  virtual void getPosMom(const StateOnPlane& state, TVector3& pos, TVector3& mom) const override;
148 
149  virtual double getMomMag(const StateOnPlane& state) const override;
150  virtual double getMomVar(const MeasuredStateOnPlane& state) const override;
151 
152  virtual TMatrixDSym get6DCov(const MeasuredStateOnPlane& state) const override;
153  virtual void getPosMomCov(const MeasuredStateOnPlane& state, TVector3& pos, TVector3& mom, TMatrixDSym& cov) const override;
154  virtual double getCharge(const StateOnPlane& state) const override;
155  virtual double getQop(const StateOnPlane& state) const override {return state.getState()(0);}
156  double getSpu(const StateOnPlane& state) const;
157  double getTime(const StateOnPlane& state) const override;
158 
159  virtual void getForwardJacobianAndNoise(TMatrixD& jacobian, TMatrixDSym& noise, TVectorD& deltaState) const override;
160 
161  virtual void getBackwardJacobianAndNoise(TMatrixD& jacobian, TMatrixDSym& noise, TVectorD& deltaState) const override;
162 
163  std::vector<genfit::MatStep> getSteps() const override;
164 
165  virtual double getRadiationLenght() const override;
166 
167  virtual void setPosMom(StateOnPlane& state, const TVector3& pos, const TVector3& mom) const override;
168  virtual void setPosMom(StateOnPlane& state, const TVectorD& state6) const override;
169  virtual void setPosMomErr(MeasuredStateOnPlane& state, const TVector3& pos, const TVector3& mom, const TVector3& posErr, const TVector3& momErr) const override;
170  virtual void setPosMomCov(MeasuredStateOnPlane& state, const TVector3& pos, const TVector3& mom, const TMatrixDSym& cov6x6) const override;
171  virtual void setPosMomCov(MeasuredStateOnPlane& state, const TVectorD& state6, const TMatrixDSym& cov6x6) const override;
172 
173  virtual void setChargeSign(StateOnPlane& state, double charge) const override;
174  virtual void setQop(StateOnPlane& state, double qop) const override {state.getState()(0) = qop;}
175 
176  void setSpu(StateOnPlane& state, double spu) const;
177  void setTime(StateOnPlane& state, double time) const override;
178 
180 
188  virtual double RKPropagate(M1x7& state7,
189  M7x7* jacobian,
190  M1x3& SA,
191  double S,
192  bool varField = true,
193  bool calcOnlyLastRowOfJ = false) const;
194 
195  virtual bool isSameType(const AbsTrackRep* other) override;
196  virtual bool isSame(const AbsTrackRep* other) override;
197 
198  private:
199 
200  void initArrays() const;
201 
202  virtual double extrapToPoint(StateOnPlane& state,
203  const TVector3& point,
204  const TMatrixDSym* G = nullptr, // weight matrix (metric)
205  bool stopAtBoundary = false,
206  bool calcJacobianNoise = false) const;
207 
208  void getState7(const StateOnPlane& state, M1x7& state7) const;
209  void getState5(StateOnPlane& state, const M1x7& state7) const; // state7 must already lie on plane of state!
210 
211  void calcJ_pM_5x7(M5x7& J_pM, const TVector3& U, const TVector3& V, const M1x3& pTilde, double spu) const;
212 
213  void transformPM6(const MeasuredStateOnPlane& state,
214  M6x6& out6x6) const;
215 
216  void calcJ_Mp_7x5(M7x5& J_Mp, const TVector3& U, const TVector3& V, const TVector3& W, const M1x3& A) const;
217 
218  void calcForwardJacobianAndNoise(const M1x7& startState7, const DetPlane& startPlane,
219  const M1x7& destState7, const DetPlane& destPlane) const;
220 
221  void transformM6P(const M6x6& in6x6,
222  const M1x7& state7,
223  MeasuredStateOnPlane& state) const; // plane and charge must already be set!
224 
226 
235  bool RKutta(const M1x4& SU,
236  const DetPlane& plane,
237  double charge,
238  double mass,
239  M1x7& state7,
240  M7x7* jacobianT,
241  M1x7* J_MMT_unprojected_lastRow,
242  double& coveredDistance, // signed
243  double& flightTime,
244  bool& checkJacProj,
245  M7x7& noiseProjection,
246  StepLimits& limits,
247  bool onlyOneStep = false,
248  bool calcOnlyLastRowOfJ = false) const;
249 
250  double estimateStep(const M1x7& state7,
251  const M1x4& SU,
252  const DetPlane& plane,
253  const double& charge,
254  double& relMomLoss,
255  StepLimits& limits) const;
256 
257  TVector3 pocaOnLine(const TVector3& linePoint,
258  const TVector3& lineDirection,
259  const TVector3& point) const;
260 
262 
270  double Extrap(const DetPlane& startPlane, // plane where Extrap starts
271  const DetPlane& destPlane, // plane where Extrap has to extrapolate to
272  double charge,
273  double mass,
274  bool& isAtBoundary,
275  M1x7& state7,
276  double& flightTime,
277  bool fillExtrapSteps,
278  TMatrixDSym* cov = nullptr,
279  bool onlyOneStep = false,
280  bool stopAtBoundary = false,
281  double maxStep = 1.E99) const;
282 
283  void checkCache(const StateOnPlane& state, const SharedPlanePtr* plane) const;
284 
285  double momMag(const M1x7& state7) const;
286 
287  protected:
288 
289  mutable StateOnPlane lastStartState_;
290  mutable StateOnPlane lastEndState_;
291 
292  private:
293 
294  mutable std::vector<RKStep> RKSteps_;
295  mutable int RKStepsFXStart_;
296  mutable int RKStepsFXStop_;
297  mutable std::vector<ExtrapStep> ExtrapSteps_;
298 
299  mutable TMatrixD fJacobian_;
300  mutable TMatrixDSym fNoise_;
301 
302  mutable bool useCache_;
303  mutable unsigned int cachePos_;
304 
305  // auxiliary variables and arrays
306  // needed in Extrap()
307  mutable StepLimits limits_;
308  mutable M7x7 noiseArray_;
309  mutable M7x7 noiseProjection_;
310  mutable M7x7 J_MMT_;
311 
312  public:
313 
314  ClassDefOverride(RKTrackRep, 1)
315 
316 };
317 
318 } /* End of namespace genfit */
321 #endif // genfit_RKTrackRep_h