Logo ROOT   6.21/01
Reference Guide
REveTrackPropagator.hxx
Go to the documentation of this file.
1 // @(#)root/eve7:$Id$
2 // Authors: Matevz Tadel & Alja Mrak-Tadel: 2006, 2007
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2019, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 #ifndef ROOT7_REveTrackPropagator
13 #define ROOT7_REveTrackPropagator
14 
15 #include <ROOT/REveVector.hxx>
16 #include <ROOT/REvePathMark.hxx>
17 #include <ROOT/REveUtil.hxx>
18 #include <ROOT/REveElement.hxx>
19 #include "TMarker.h"
20 
21 #include <vector>
22 
23 namespace ROOT {
24 namespace Experimental {
25 
26 class REvePointSet;
27 
28 ////////////////////////////////////////////////////////////////////////////////
29 /// REveMagField
30 /// Abstract interface to magnetic field
31 ////////////////////////////////////////////////////////////////////////////////
32 
34 {
35 protected:
37 
38 public:
39  REveMagField() = default;
40  virtual ~REveMagField() {}
41 
42  virtual Bool_t IsConst() const { return fFieldConstant; }
43 
44  virtual void PrintField(Double_t x, Double_t y, Double_t z) const
45  {
46  REveVector b = GetField(x, y, z);
47  printf("v(%f, %f, %f) B(%f, %f, %f) \n", x, y, z, b.fX, b.fY, b.fZ);
48  }
49 
50  REveVectorD GetFieldD(const REveVectorD &v) const { return GetFieldD(v.fX, v.fY, v.fZ); }
51 
52  // Track propgator uses only GetFieldD() and GetMaxFieldMagD(). Have to keep/reuse
53  // GetField() and GetMaxFieldMag() because of backward compatibility.
54 
55  virtual REveVectorD GetFieldD(Double_t x, Double_t y, Double_t z) const { return GetField(x, y, z); }
56  virtual Double_t GetMaxFieldMagD() const
57  {
58  return GetMaxFieldMag();
59  } // not abstract because of backward compatibility
60 
61  virtual REveVector GetField(Float_t, Float_t, Float_t) const { return REveVector(); }
62  virtual Float_t GetMaxFieldMag() const { return 4; } // not abstract because of backward compatibility
63 
64 };
65 
66 ////////////////////////////////////////////////////////////////////////////////
67 /// REveMagFieldConst
68 /// Interface to constant magnetic field.
69 ////////////////////////////////////////////////////////////////////////////////
70 
72 {
73 protected:
75 
76 public:
78  virtual ~REveMagFieldConst() {}
79 
80  REveVectorD GetFieldD(Double_t /*x*/, Double_t /*y*/, Double_t /*z*/) const override { return fB; }
81 
82  Double_t GetMaxFieldMagD() const override { return fB.Mag(); };
83 };
84 
85 ////////////////////////////////////////////////////////////////////////////////
86 /// REveMagFieldDuo
87 /// Interface to magnetic field with two different values depending on radius.
88 ////////////////////////////////////////////////////////////////////////////////
89 
91 {
92 protected:
96 
97 public:
99  : REveMagField(), fBIn(0, 0, bIn), fBOut(0, 0, bOut), fR2(r * r)
100  {
102  }
103  virtual ~REveMagFieldDuo() {}
104 
106  {
107  return ((x * x + y * y) < fR2) ? fBIn : fBOut;
108  }
109 
110  Double_t GetMaxFieldMagD() const override
111  {
112  Double_t b1 = fBIn.Mag(), b2 = fBOut.Mag();
113  return b1 > b2 ? b1 : b2;
114  }
115 };
116 
117 ////////////////////////////////////////////////////////////////////////////////
118 /// REveTrackPropagator
119 /// Calculates path of a particle taking into account special path-marks and imposed boundaries.
120 ////////////////////////////////////////////////////////////////////////////////
121 
123  public REveRefBackPtr
124 {
125 public:
127 
129 
130 protected:
131  struct Helix_t {
132  Int_t fCharge; // Charge of tracked particle.
133  Double_t fMaxAng; // Maximum step angle.
134  Double_t fMaxStep; // Maximum allowed step size.
135  Double_t fDelta; // Maximum error in the middle of the step.
136 
137  Double_t fPhi; // Accumulated angle to check fMaxOrbs by propagator.
138  Bool_t fValid; // Corner case pT~0 or B~0, possible in variable mag field.
139 
140  // ----------------------------------------------------------------
141 
142  // helix parameters
143  Double_t fLam; // Momentum ratio pT/pZ.
144  Double_t fR; // Helix radius in cm.
145  Double_t fPhiStep; // Caluclated from fMinAng and fDelta.
146  Double_t fSin, fCos; // Current sin/cos(phistep).
147 
148  // Runge-Kutta parameters
149  Double_t fRKStep; // Step for Runge-Kutta.
150 
151  // cached
152  REveVectorD fB; // Current magnetic field, cached.
153  REveVectorD fE1, fE2, fE3; // Base vectors: E1 -> B dir, E2->pT dir, E3 = E1xE2.
154  REveVectorD fPt, fPl; // Transverse and longitudinal momentum.
155  Double_t fPtMag; // Magnitude of pT.
156  Double_t fPlMag; // Momentum parallel to mag field.
157  Double_t fLStep; // Transverse step arc-length in cm.
158 
159  // ----------------------------------------------------------------
160 
161  Helix_t();
162 
163  void UpdateCommon(const REveVectorD &p, const REveVectorD &b);
164  void UpdateHelix(const REveVectorD &p, const REveVectorD &b, Bool_t full_update, Bool_t enforce_max_step);
165  void UpdateRK(const REveVectorD &p, const REveVectorD &b);
166 
167  void Step(const REveVector4D &v, const REveVectorD &p, REveVector4D &vOut, REveVectorD &pOut);
168 
169  Double_t GetStep() { return fLStep * TMath::Sqrt(1 + fLam * fLam); }
170  Double_t GetStep2() { return fLStep * fLStep * (1 + fLam * fLam); }
171  };
172 
173 private:
174  REveTrackPropagator(const REveTrackPropagator &) = delete;
176 
177  void DistributeOffset(const REveVectorD &off, Int_t first_point, Int_t np, REveVectorD &p);
178 
179 protected:
181 
184 
185  // Track extrapolation limits
186  Double_t fMaxR; // Max radius for track extrapolation
187  Double_t fMaxZ; // Max z-coordinate for track extrapolation.
188  Int_t fNMax; // Max steps
189  // Helix limits
190  Double_t fMaxOrbs; // Maximal angular path of tracks' orbits (1 ~ 2Pi).
191 
192  // Path-mark / first-vertex control
193  Bool_t fEditPathMarks; // Show widgets for path-mark control in GUI editor.
194  Bool_t fFitDaughters; // Pass through daughter creation points when extrapolating a track.
195  Bool_t fFitReferences; // Pass through given track-references when extrapolating a track.
196  Bool_t fFitDecay; // Pass through decay point when extrapolating a track.
197  Bool_t fFitCluster2Ds; // Pass through 2D-clusters when extrapolating a track.
198  Bool_t fFitLineSegments; // Pass through line when extrapolating a track.
199  Bool_t fRnrDaughters; // Render daughter path-marks.
200  Bool_t fRnrReferences; // Render track-reference path-marks.
201  Bool_t fRnrDecay; // Render decay path-marks.
202  Bool_t fRnrCluster2Ds; // Render 2D-clusters.
203  Bool_t fRnrFV; // Render first vertex.
204  TMarker fPMAtt; // Marker attributes for rendering of path-marks.
205  TMarker fFVAtt; // Marker attributes for fits vertex.
206 
207  // Handling of discontinuities in projections
208  UChar_t fProjTrackBreaking; // Handling of projected-track breaking.
209  Bool_t fRnrPTBMarkers; // Render break-points on tracks.
210  TMarker fPTBAtt; // Marker attributes for track break-points.
211 
212  // ----------------------------------------------------------------
213 
214  // Propagation, state of current track
215  std::vector<REveVector4D> fPoints; // Calculated point.
216  std::vector<REveVector4D> fLastPoints; // Copy of the latest calculated points.
217  REveVectorD fV; // Start vertex.
218  Helix_t fH; // Helix.
219 
220  void RebuildTracks();
221  void Update(const REveVector4D &v, const REveVectorD &p, Bool_t full_update = kFALSE, Bool_t enforce_max_step = kFALSE);
222  void Step(const REveVector4D &v, const REveVectorD &p, REveVector4D &vOut, REveVectorD &pOut);
223 
226  void LoopToBounds(REveVectorD &p);
227 
229  void LineToBounds(REveVectorD &p);
230 
231  void StepRungeKutta(Double_t step, Double_t *vect, Double_t *vout);
232 
233  Bool_t HelixIntersectPlane(const REveVectorD &p, const REveVectorD &point, const REveVectorD &normal, REveVectorD &itsect);
234  Bool_t LineIntersectPlane(const REveVectorD &p, const REveVectorD &point, const REveVectorD &normal, REveVectorD &itsect);
235  Bool_t PointOverVertex(const REveVector4D &v0, const REveVector4D &v, Double_t *p = 0);
236 
238  Double_t rMagInv, REveVectorD &c);
240  REveVectorD &out);
241 
242 public:
243  REveTrackPropagator(const std::string& n = "REveTrackPropagator", const std::string& t = "", REveMagField *field = nullptr,
244  Bool_t own_field = kTRUE);
245  virtual ~REveTrackPropagator();
246 
247  void OnZeroRefCount() override;
248 
249  void CheckReferenceCount(const std::string &from = "<unknown>") override;
250 
251  void StampAllTracks();
252 
253  // propagation
254  void InitTrack(const REveVectorD &v, Int_t charge);
255  void ResetTrack();
256 
257  Int_t GetCurrentPoint() const;
258  Double_t GetTrackLength(Int_t start_point = 0, Int_t end_point = -1) const;
259 
260  virtual void GoToBounds(REveVectorD &p);
261  virtual Bool_t GoToVertex(REveVectorD &v, REveVectorD &p);
262  virtual Bool_t GoToLineSegment(const REveVectorD &s, const REveVectorD &r, REveVectorD &p);
263 
264  // REveVectorF wrappers
265  void InitTrack(const REveVectorF &v, Int_t charge);
266  void GoToBounds(REveVectorF &p);
269 
270  Bool_t IntersectPlane(const REveVectorD &p, const REveVectorD &point, const REveVectorD &normal, REveVectorD &itsect);
271 
272  void FillPointSet(REvePointSet *ps) const;
273 
275 
276  void SetMagField(Double_t bX, Double_t bY, Double_t bZ);
277  void SetMagField(Double_t b) { SetMagField(0, 0, b); }
278  void SetMagFieldObj(REveMagField *field, Bool_t own_field = kTRUE);
279 
280  void SetMaxR(Double_t x);
281  void SetMaxZ(Double_t x);
282  void SetMaxOrbs(Double_t x);
283  void SetMinAng(Double_t x);
284  void SetMaxAng(Double_t x);
285  void SetMaxStep(Double_t x);
286  void SetDelta(Double_t x);
287 
289  void SetRnrDaughters(Bool_t x);
290  void SetRnrReferences(Bool_t x);
291  void SetRnrDecay(Bool_t x);
292  void SetRnrCluster2Ds(Bool_t x);
293  void SetFitDaughters(Bool_t x);
294  void SetFitReferences(Bool_t x);
295  void SetFitDecay(Bool_t x);
296  void SetFitCluster2Ds(Bool_t x);
298  void SetRnrFV(Bool_t x);
300  void SetRnrPTBMarkers(Bool_t x);
301 
303  void PrintMagField(Double_t x, Double_t y, Double_t z) const;
304 
305  EStepper_e GetStepper() const { return fStepper; }
306 
307  Double_t GetMaxR() const { return fMaxR; }
308  Double_t GetMaxZ() const { return fMaxZ; }
309  Double_t GetMaxOrbs() const { return fMaxOrbs; }
310  Double_t GetMinAng() const;
311  Double_t GetMaxAng() const { return fH.fMaxAng; }
312  Double_t GetMaxStep() const { return fH.fMaxStep; }
313  Double_t GetDelta() const { return fH.fDelta; }
314 
318  Bool_t GetRnrDecay() const { return fRnrDecay; }
322  Bool_t GetFitDecay() const { return fFitDecay; }
325  Bool_t GetRnrFV() const { return fRnrFV; }
328 
329  TMarker &RefPMAtt() { return fPMAtt; }
330  TMarker &RefFVAtt() { return fFVAtt; }
331  TMarker &RefPTBAtt() { return fPTBAtt; }
332 
333  const std::vector<REveVector4D> &GetLastPoints() const { return fLastPoints; }
334 
335  static Bool_t IsOutsideBounds(const REveVectorD &point, Double_t maxRsqr, Double_t maxZ);
336 
337  static Double_t fgDefMagField; // Default value for constant solenoid magnetic field.
338  static const Double_t fgkB2C; // Constant for conversion of momentum to curvature.
339  static REveTrackPropagator fgDefault; // Default track propagator.
340 
341  static Double_t fgEditorMaxR; // Max R that can be set in GUI editor.
342  static Double_t fgEditorMaxZ; // Max Z that can be set in GUI editor.
343 };
344 
345 //______________________________________________________________________________
347 {
348  // Return true if point% is outside of cylindrical bounds detrmined by
349  // square radius and z.
350 
351  return TMath::Abs(point.fZ) > maxZ || point.fX * point.fX + point.fY * point.fY > maxRsqr;
352 }
353 
354 //______________________________________________________________________________
356 {
357  static const Double_t kMinPl = 1e-5;
358 
359  REveVectorD dv;
360  dv.Sub(v0, v);
361 
362  Double_t dotV;
363 
364  if (TMath::Abs(fH.fPlMag) > kMinPl) {
365  // Use longitudinal momentum to determine crossing point.
366  // Works ok for spiraling helices, also for loopers.
367 
368  dotV = fH.fE1.Dot(dv);
369  if (fH.fPlMag < 0)
370  dotV = -dotV;
371  } else {
372  // Use full momentum, which is pT, under this conditions.
373 
374  dotV = fH.fE2.Dot(dv);
375  }
376 
377  if (p)
378  *p = dotV;
379 
380  return dotV < 0;
381 }
382 
383 } // namespace Experimental
384 } // namespace ROOT
385 
386 #endif
void SetProjTrackBreaking(UChar_t x)
Set projection break-point mode and rebuild tracks.
Bool_t HelixIntersectPlane(const REveVectorD &p, const REveVectorD &point, const REveVectorD &normal, REveVectorD &itsect)
Intersect helix with a plane.
virtual void GoToBounds(REveVectorD &p)
Propagate particle to bounds.
REveVectorD GetFieldD(const REveVectorD &v) const
void InitTrack(const REveVectorD &v, Int_t charge)
Initialize internal data-members for given particle parameters.
void CheckReferenceCount(const std::string &from="<unknown>") override
Check reference count - virtual from REveElement.
REveTrackPropagator Calculates path of a particle taking into account special path-marks and imposed ...
void DistributeOffset(const REveVectorD &off, Int_t first_point, Int_t np, REveVectorD &p)
Distribute offset between first and last point index and rotate momentum.
Bool_t IntersectPlane(const REveVectorD &p, const REveVectorD &point, const REveVectorD &normal, REveVectorD &itsect)
Find intersection of currently propagated track with a plane.
Returns the available number of logical cores.
Definition: StringConv.hxx:21
float Float_t
Definition: RtypesCore.h:53
void RebuildTracks()
Rebuild all tracks using this render-style.
virtual Double_t GetMaxFieldMagD() const
virtual REveVectorD GetFieldD(Double_t x, Double_t y, Double_t z) const
const std::vector< REveVector4D > & GetLastPoints() const
Bool_t ClosestPointBetweenLines(const REveVectorD &, const REveVectorD &, const REveVectorD &, const REveVectorD &, REveVectorD &out)
Get closest point on line defined with vector p0 and u.
Double_t GetMinAng() const
Get maximum step angle.
virtual void PrintField(Double_t x, Double_t y, Double_t z) const
REveVectorD GetFieldD(Double_t, Double_t, Double_t) const override
REveTrackPropagator & operator=(const REveTrackPropagator &)=delete
void SetRnrReferences(Bool_t x)
Set track-reference rendering and rebuild tracks.
static constexpr double ps
Manages Markers.
Definition: TMarker.h:23
REveRefBackPtr reference-count with back pointers.
Definition: REveUtil.hxx:130
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
void SetRnrFV(Bool_t x)
Set first-vertex rendering and rebuild tracks.
void SetDelta(Double_t x)
Set maximum error and rebuild tracks.
Int_t GetCurrentPoint() const
Get index of current point on track.
void FillPointSet(REvePointSet *ps) const
Reset ps and populate it with points in propagation cache.
void UpdateCommon(const REveVectorD &p, const REveVectorD &b)
Common update code for helix and RK propagation.
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
void SetFitReferences(Bool_t x)
Set track-reference fitting and rebuild tracks.
void LoopToBounds(REveVectorD &p)
Propagate charged particle with momentum p to bounds.
void Update(const REveVector4D &v, const REveVectorD &p, Bool_t full_update=kFALSE, Bool_t enforce_max_step=kFALSE)
Update helix / B-field projection state.
Double_t x[n]
Definition: legend1.C:17
void SetMinAng(Double_t x)
Set maximum step angle and rebuild tracks.
virtual REveVector GetField(Float_t, Float_t, Float_t) const
void SetFitDaughters(Bool_t x)
Set daughter creation point fitting and rebuild tracks.
void SetRnrDaughters(Bool_t x)
Set daughter rendering and rebuild tracks.
static constexpr double s
REveVectorD GetFieldD(Double_t x, Double_t y, Double_t) const override
void SetMagField(Double_t bX, Double_t bY, Double_t bZ)
Set constant magnetic field and rebuild tracks.
REveVectorT< Float_t > REveVector
Definition: REveVector.hxx:119
Bool_t LoopToVertex(REveVectorD &v, REveVectorD &p)
Propagate charged particle with momentum p to vertex v.
virtual Bool_t GoToVertex(REveVectorD &v, REveVectorD &p)
Propagate particle with momentum p to vertex v.
TT Dot(const REveVectorT &a) const
Definition: REveVector.hxx:164
void LineToBounds(REveVectorD &p)
Propagate neutral particle with momentum p to bounds.
REveMagFieldDuo Interface to magnetic field with two different values depending on radius...
REveVectorD GetMagField(Double_t x, Double_t y, Double_t z)
REveMagField Abstract interface to magnetic field.
ROOT::R::TRInterface & r
Definition: Object.C:4
void SetMaxZ(Double_t x)
Set maximum z and rebuild tracks.
Bool_t PointOverVertex(const REveVector4D &v0, const REveVector4D &v, Double_t *p=0)
void SetMaxStep(Double_t x)
Set maximum step-size and rebuild tracks.
void ClosestPointFromVertexToLineSegment(const REveVectorD &v, const REveVectorD &s, const REveVectorD &r, Double_t rMagInv, REveVectorD &c)
Get closest point from given vertex v to line segment defined with s and r.
void ResetTrack()
Reset cache holding particle trajectory.
void StepRungeKutta(Double_t step, Double_t *vect, Double_t *vout)
Wrapper to step with method RungeKutta.
void SetMaxOrbs(Double_t x)
Set maximum number of orbits and rebuild tracks.
void SetRnrDecay(Bool_t x)
Set decay rendering and rebuild tracks.
Bool_t LineIntersectPlane(const REveVectorD &p, const REveVectorD &point, const REveVectorD &normal, REveVectorD &itsect)
Intersect line with a plane.
void UpdateHelix(const REveVectorD &p, const REveVectorD &b, Bool_t full_update, Bool_t enforce_max_step)
Update helix parameters.
virtual Bool_t GoToLineSegment(const REveVectorD &s, const REveVectorD &r, REveVectorD &p)
Propagate particle with momentum p to line with start point s and vector r to the second point...
const Bool_t kFALSE
Definition: RtypesCore.h:88
static Bool_t IsOutsideBounds(const REveVectorD &point, Double_t maxRsqr, Double_t maxZ)
void SetFitLineSegments(Bool_t x)
Set line segment fitting and rebuild tracks.
void SetFitDecay(Bool_t x)
Set decay fitting and rebuild tracks.
void SetMaxAng(Double_t x)
Set maximum step angle and rebuild tracks.
void SetMagFieldObj(REveMagField *field, Bool_t own_field=kTRUE)
Set constant magnetic field and rebuild tracks.
Bool_t LoopToLineSegment(const REveVectorD &s, const REveVectorD &r, REveVectorD &p)
Propagate charged particle with momentum p to line segment with point s and vector r to the second po...
double Double_t
Definition: RtypesCore.h:55
void Step(const REveVector4D &v, const REveVectorD &p, REveVector4D &vOut, REveVectorD &pOut)
Step helix for given momentum p from vertex v.
void SetRnrPTBMarkers(Bool_t x)
Set projection break-point rendering and rebuild tracks.
Double_t y[n]
Definition: legend1.C:17
Bool_t LineToVertex(REveVectorD &v)
Propagate neutral particle to vertex v.
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
void SetRnrCluster2Ds(Bool_t x)
Set rendering of 2D-clusters and rebuild tracks.
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
Double_t GetTrackLength(Int_t start_point=0, Int_t end_point=-1) const
Calculate track length from start_point to end_point.
void PrintMagField(Double_t x, Double_t y, Double_t z) const
void Step(const REveVector4D &v, const REveVectorD &p, REveVector4D &vOut, REveVectorD &pOut)
Wrapper to step helix.
REveTrackPropagator(const REveTrackPropagator &)=delete
Double_t GetMaxFieldMagD() const override
virtual Float_t GetMaxFieldMag() const
void SetFitCluster2Ds(Bool_t x)
Set 2D-cluster fitting and rebuild tracks.
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
void StampAllTracks()
Element-change notification.
#define c(i)
Definition: RSha256.hxx:101
REveMagFieldConst Interface to constant magnetic field.
unsigned char UChar_t
Definition: RtypesCore.h:34
Double_t Sqrt(Double_t x)
Definition: TMath.h:681
const Bool_t kTRUE
Definition: RtypesCore.h:87
REveVector4T A four-vector template without TObject inheritance and virtual functions.
Definition: REveVector.hxx:238
const Int_t n
Definition: legend1.C:16
void UpdateRK(const REveVectorD &p, const REveVectorD &b)
Update helix for stepper RungeKutta.
void OnZeroRefCount() override
Virtual from REveRefBackPtr - track reference count has reached zero.
REveVectorT & Sub(const REveVectorT &a, const REveVectorT &b)
Definition: REveVector.hxx:182
void SetMaxR(Double_t x)
Set maximum radius and rebuild tracks.
REveMagFieldDuo(Double_t r, Double_t bIn, Double_t bOut)
REveMagFieldConst(Double_t x, Double_t y, Double_t z)