Logo ROOT  
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
23namespace ROOT {
24namespace Experimental {
25
26class REvePointSet;
27
28////////////////////////////////////////////////////////////////////////////////
29/// REveMagField
30/// Abstract interface to magnetic field
31////////////////////////////////////////////////////////////////////////////////
32
34{
35protected:
37
38public:
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 virtual Double_t GetMaxFieldMag() const = 0;
52
53 REveVectorD GetField(const REveVectorD &v) const { return GetField(v.fX, v.fY, v.fZ); }
54};
55
56////////////////////////////////////////////////////////////////////////////////
57/// REveMagFieldConst
58/// Interface to constant magnetic field.
59////////////////////////////////////////////////////////////////////////////////
60
62{
63protected:
65
66public:
68 virtual ~REveMagFieldConst() {}
69
70 Double_t GetMaxFieldMag() const override { return fB.Mag(); };
71 REveVectorD GetField(Double_t /*x*/, Double_t /*y*/, Double_t /*z*/) const override { return fB; }
72};
73
74////////////////////////////////////////////////////////////////////////////////
75/// REveMagFieldDuo
76/// Interface to magnetic field with two different values depending on radius.
77////////////////////////////////////////////////////////////////////////////////
78
80{
81protected:
85
86public:
88 : REveMagField(), fBIn(0, 0, bIn), fBOut(0, 0, bOut), fR2(r * r)
89 {
91 }
92 virtual ~REveMagFieldDuo() {}
93
94 Double_t GetMaxFieldMag() const override
95 {
96 Double_t b1 = fBIn.Mag(), b2 = fBOut.Mag();
97 return b1 > b2 ? b1 : b2;
98 }
99
101 {
102 return ((x * x + y * y) < fR2) ? fBIn : fBOut;
103 }
104};
105
106////////////////////////////////////////////////////////////////////////////////
107/// REveTrackPropagator
108/// Calculates path of a particle taking into account special path-marks and imposed boundaries.
109////////////////////////////////////////////////////////////////////////////////
110
112 public REveRefBackPtr
113{
114public:
116
118
119protected:
120 struct Helix_t {
121 Int_t fCharge; // Charge of tracked particle.
122 Double_t fMaxAng; // Maximum step angle.
123 Double_t fMaxStep; // Maximum allowed step size.
124 Double_t fDelta; // Maximum error in the middle of the step.
125
126 Double_t fPhi; // Accumulated angle to check fMaxOrbs by propagator.
127 Bool_t fValid; // Corner case pT~0 or B~0, possible in variable mag field.
128
129 // ----------------------------------------------------------------
130
131 // helix parameters
132 Double_t fLam; // Momentum ratio pT/pZ.
133 Double_t fR; // Helix radius in cm.
134 Double_t fPhiStep; // Caluclated from fMinAng and fDelta.
135 Double_t fSin, fCos; // Current sin/cos(phistep).
136
137 // Runge-Kutta parameters
138 Double_t fRKStep; // Step for Runge-Kutta.
139
140 // cached
141 REveVectorD fB; // Current magnetic field, cached.
142 REveVectorD fE1, fE2, fE3; // Base vectors: E1 -> B dir, E2->pT dir, E3 = E1xE2.
143 REveVectorD fPt, fPl; // Transverse and longitudinal momentum.
144 Double_t fPtMag; // Magnitude of pT.
145 Double_t fPlMag; // Momentum parallel to mag field.
146 Double_t fLStep; // Transverse step arc-length in cm.
147
148 // ----------------------------------------------------------------
149
150 Helix_t();
151
152 void UpdateCommon(const REveVectorD &p, const REveVectorD &b);
153 void UpdateHelix(const REveVectorD &p, const REveVectorD &b, Bool_t full_update, Bool_t enforce_max_step);
154 void UpdateRK(const REveVectorD &p, const REveVectorD &b);
155
156 void Step(const REveVector4D &v, const REveVectorD &p, REveVector4D &vOut, REveVectorD &pOut);
157
159 Double_t GetStep2() { return fLStep * fLStep * (1 + fLam * fLam); }
160 };
161
162private:
165
166 void DistributeOffset(const REveVectorD &off, Int_t first_point, Int_t np, REveVectorD &p);
167
168protected:
170
173
174 // Track extrapolation limits
175 Double_t fMaxR; // Max radius for track extrapolation
176 Double_t fMaxZ; // Max z-coordinate for track extrapolation.
177 Int_t fNMax; // Max steps
178 // Helix limits
179 Double_t fMaxOrbs; // Maximal angular path of tracks' orbits (1 ~ 2Pi).
180
181 // Path-mark / first-vertex control
182 Bool_t fEditPathMarks; // Show widgets for path-mark control in GUI editor.
183 Bool_t fFitDaughters; // Pass through daughter creation points when extrapolating a track.
184 Bool_t fFitReferences; // Pass through given track-references when extrapolating a track.
185 Bool_t fFitDecay; // Pass through decay point when extrapolating a track.
186 Bool_t fFitCluster2Ds; // Pass through 2D-clusters when extrapolating a track.
187 Bool_t fFitLineSegments; // Pass through line when extrapolating a track.
188 Bool_t fRnrDaughters; // Render daughter path-marks.
189 Bool_t fRnrReferences; // Render track-reference path-marks.
190 Bool_t fRnrDecay; // Render decay path-marks.
191 Bool_t fRnrCluster2Ds; // Render 2D-clusters.
192 Bool_t fRnrFV; // Render first vertex.
193 TMarker fPMAtt; // Marker attributes for rendering of path-marks.
194 TMarker fFVAtt; // Marker attributes for fits vertex.
195
196 // Handling of discontinuities in projections
197 UChar_t fProjTrackBreaking; // Handling of projected-track breaking.
198 Bool_t fRnrPTBMarkers; // Render break-points on tracks.
199 TMarker fPTBAtt; // Marker attributes for track break-points.
200
201 // ----------------------------------------------------------------
202
203 // Propagation, state of current track
204 std::vector<REveVector4D> fPoints; // Calculated point.
205 std::vector<REveVector4D> fLastPoints; // Copy of the latest calculated points.
206 REveVectorD fV; // Start vertex.
207 Helix_t fH; // Helix.
208
209 void RebuildTracks();
210 void Update(const REveVector4D &v, const REveVectorD &p, Bool_t full_update = kFALSE, Bool_t enforce_max_step = kFALSE);
211 void Step(const REveVector4D &v, const REveVectorD &p, REveVector4D &vOut, REveVectorD &pOut);
212
215 void LoopToBounds(REveVectorD &p);
216
218 void LineToBounds(REveVectorD &p);
219
220 void StepRungeKutta(Double_t step, Double_t *vect, Double_t *vout);
221
222 Bool_t HelixIntersectPlane(const REveVectorD &p, const REveVectorD &point, const REveVectorD &normal, REveVectorD &itsect);
223 Bool_t LineIntersectPlane(const REveVectorD &p, const REveVectorD &point, const REveVectorD &normal, REveVectorD &itsect);
225
227 Double_t rMagInv, REveVectorD &c);
229 REveVectorD &out);
230
231public:
232 REveTrackPropagator(const std::string& n = "REveTrackPropagator", const std::string& t = "", REveMagField *field = nullptr,
233 Bool_t own_field = kTRUE);
234 virtual ~REveTrackPropagator();
235
236 void OnZeroRefCount() override;
237
238 void CheckReferenceCount(const std::string &from = "<unknown>") override;
239
240 void StampAllTracks();
241
242 // propagation
243 void InitTrack(const REveVectorD &v, Int_t charge);
244 void ResetTrack();
245
246 Int_t GetCurrentPoint() const;
247 Double_t GetTrackLength(Int_t start_point = 0, Int_t end_point = -1) const;
248
249 virtual void GoToBounds(REveVectorD &p);
251 virtual Bool_t GoToLineSegment(const REveVectorD &s, const REveVectorD &r, REveVectorD &p);
252
253 // REveVectorF wrappers
254 void InitTrack(const REveVectorF &v, Int_t charge);
255 void GoToBounds(REveVectorF &p);
258
259 Bool_t IntersectPlane(const REveVectorD &p, const REveVectorD &point, const REveVectorD &normal, REveVectorD &itsect);
260
261 void FillPointSet(REvePointSet *ps) const;
262
264
265 void SetMagField(Double_t bX, Double_t bY, Double_t bZ);
267 void SetMagFieldObj(REveMagField *field, Bool_t own_field = kTRUE);
268
269 void SetMaxR(Double_t x);
270 void SetMaxZ(Double_t x);
271 void SetMaxOrbs(Double_t x);
272 void SetMinAng(Double_t x);
273 void SetMaxAng(Double_t x);
274 void SetMaxStep(Double_t x);
275 void SetDelta(Double_t x);
276
280 void SetRnrDecay(Bool_t x);
284 void SetFitDecay(Bool_t x);
287 void SetRnrFV(Bool_t x);
290
292 void PrintMagField(Double_t x, Double_t y, Double_t z) const;
293
294 EStepper_e GetStepper() const { return fStepper; }
295
296 Double_t GetMaxR() const { return fMaxR; }
297 Double_t GetMaxZ() const { return fMaxZ; }
298 Double_t GetMaxOrbs() const { return fMaxOrbs; }
299 Double_t GetMinAng() const;
300 Double_t GetMaxAng() const { return fH.fMaxAng; }
301 Double_t GetMaxStep() const { return fH.fMaxStep; }
302 Double_t GetDelta() const { return fH.fDelta; }
303
307 Bool_t GetRnrDecay() const { return fRnrDecay; }
311 Bool_t GetFitDecay() const { return fFitDecay; }
314 Bool_t GetRnrFV() const { return fRnrFV; }
317
318 TMarker &RefPMAtt() { return fPMAtt; }
319 TMarker &RefFVAtt() { return fFVAtt; }
320 TMarker &RefPTBAtt() { return fPTBAtt; }
321
322 const std::vector<REveVector4D> &GetLastPoints() const { return fLastPoints; }
323
324 static Bool_t IsOutsideBounds(const REveVectorD &point, Double_t maxRsqr, Double_t maxZ);
325
326 static Double_t fgDefMagField; // Default value for constant solenoid magnetic field.
327 static const Double_t fgkB2C; // Constant for conversion of momentum to curvature.
328 static REveTrackPropagator fgDefault; // Default track propagator.
329};
330
331//______________________________________________________________________________
333{
334 // Return true if point% is outside of cylindrical bounds detrmined by
335 // square radius and z.
336
337 return TMath::Abs(point.fZ) > maxZ || point.fX * point.fX + point.fY * point.fY > maxRsqr;
338}
339
340//______________________________________________________________________________
342{
343 static const Double_t kMinPl = 1e-5;
344
345 REveVectorD dv;
346 dv.Sub(v0, v);
347
348 Double_t dotV;
349
350 if (TMath::Abs(fH.fPlMag) > kMinPl) {
351 // Use longitudinal momentum to determine crossing point.
352 // Works ok for spiraling helices, also for loopers.
353
354 dotV = fH.fE1.Dot(dv);
355 if (fH.fPlMag < 0)
356 dotV = -dotV;
357 } else {
358 // Use full momentum, which is pT, under this conditions.
359
360 dotV = fH.fE2.Dot(dv);
361 }
362
363 if (p)
364 *p = dotV;
365
366 return dotV < 0;
367}
368
369} // namespace Experimental
370} // namespace ROOT
371
372#endif
ROOT::R::TRInterface & r
Definition: Object.C:4
#define b(i)
Definition: RSha256.hxx:100
#define c(i)
Definition: RSha256.hxx:101
#define e(i)
Definition: RSha256.hxx:103
int Int_t
Definition: RtypesCore.h:45
unsigned char UChar_t
Definition: RtypesCore.h:38
const Bool_t kFALSE
Definition: RtypesCore.h:101
bool Bool_t
Definition: RtypesCore.h:63
double Double_t
Definition: RtypesCore.h:59
const Bool_t kTRUE
Definition: RtypesCore.h:100
REvePointSet is a render-element holding a collection of 3D points with optional per-point TRef and a...
REveMagFieldConst Interface to constant magnetic field.
REveVectorD GetField(Double_t, Double_t, Double_t) const override
REveMagFieldConst(Double_t x, Double_t y, Double_t z)
REveMagFieldDuo Interface to magnetic field with two different values depending on radius.
REveVectorD GetField(Double_t x, Double_t y, Double_t) const override
REveMagFieldDuo(Double_t r, Double_t bIn, Double_t bOut)
Double_t GetMaxFieldMag() const override
REveMagField Abstract interface to magnetic field.
virtual Double_t GetMaxFieldMag() const =0
REveVectorD GetField(const REveVectorD &v) const
virtual void PrintField(Double_t x, Double_t y, Double_t z) const
virtual REveVectorD GetField(Double_t x, Double_t y, Double_t z) const =0
REveRefBackPtr reference-count with back pointers.
Definition: REveUtil.hxx:130
REveTrackPropagator Calculates path of a particle taking into account special path-marks and imposed ...
virtual Bool_t GoToVertex(REveVectorD &v, REveVectorD &p)
Propagate particle with momentum p to vertex v.
REveTrackPropagator & operator=(const REveTrackPropagator &)=delete
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.
void SetRnrPTBMarkers(Bool_t x)
Set projection break-point rendering and rebuild tracks.
REveVectorD GetMagField(Double_t x, Double_t y, Double_t z)
static Bool_t IsOutsideBounds(const REveVectorD &point, Double_t maxRsqr, Double_t maxZ)
void SetRnrDaughters(Bool_t x)
Set daughter rendering and rebuild tracks.
Bool_t HelixIntersectPlane(const REveVectorD &p, const REveVectorD &point, const REveVectorD &normal, REveVectorD &itsect)
Intersect helix with a plane.
void SetFitLineSegments(Bool_t x)
Set line segment fitting and rebuild tracks.
void SetRnrCluster2Ds(Bool_t x)
Set rendering of 2D-clusters and rebuild tracks.
Bool_t LineIntersectPlane(const REveVectorD &p, const REveVectorD &point, const REveVectorD &normal, REveVectorD &itsect)
Intersect line with a plane.
void CheckReferenceCount(const std::string &from="<unknown>") override
Check reference count - virtual from REveElement.
Double_t GetMinAng() const
Get maximum step angle.
void ResetTrack()
Reset cache holding particle trajectory.
void Step(const REveVector4D &v, const REveVectorD &p, REveVector4D &vOut, REveVectorD &pOut)
Wrapper to step helix.
Double_t GetTrackLength(Int_t start_point=0, Int_t end_point=-1) const
Calculate track length from start_point to end_point.
void InitTrack(const REveVectorD &v, Int_t charge)
Initialize internal data-members for given particle parameters.
void LineToBounds(REveVectorD &p)
Propagate neutral particle with momentum p to bounds.
void SetRnrReferences(Bool_t x)
Set track-reference rendering 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.
Bool_t IntersectPlane(const REveVectorD &p, const REveVectorD &point, const REveVectorD &normal, REveVectorD &itsect)
Find intersection of currently propagated track with a plane.
void RebuildTracks()
Rebuild all tracks using this render-style.
void SetMaxAng(Double_t x)
Set maximum step angle and rebuild tracks.
Int_t GetCurrentPoint() const
Get index of current point on track.
void SetMaxZ(Double_t x)
Set maximum z and rebuild tracks.
void SetMaxStep(Double_t x)
Set maximum step-size and rebuild tracks.
void OnZeroRefCount() override
Virtual from REveRefBackPtr - track reference count has reached zero.
void SetRnrDecay(Bool_t x)
Set decay rendering and rebuild tracks.
Bool_t PointOverVertex(const REveVector4D &v0, const REveVector4D &v, Double_t *p=0)
void FillPointSet(REvePointSet *ps) const
Reset ps and populate it with points in propagation cache.
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...
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.
const std::vector< REveVector4D > & GetLastPoints() const
void SetFitDaughters(Bool_t x)
Set daughter creation point fitting and rebuild tracks.
void SetFitDecay(Bool_t x)
Set decay fitting and rebuild tracks.
void SetFitReferences(Bool_t x)
Set track-reference fitting and rebuild tracks.
void StampAllTracks()
Element-change notification.
REveTrackPropagator(const REveTrackPropagator &)=delete
void SetMaxOrbs(Double_t x)
Set maximum number of orbits and rebuild tracks.
void StepRungeKutta(Double_t step, Double_t *vect, Double_t *vout)
Wrapper to step with method RungeKutta.
void SetProjTrackBreaking(UChar_t x)
Set projection break-point mode and rebuild tracks.
Bool_t ClosestPointBetweenLines(const REveVectorD &, const REveVectorD &, const REveVectorD &, const REveVectorD &, REveVectorD &out)
Get closest point on line defined with vector p0 and u.
void SetMagFieldObj(REveMagField *field, Bool_t own_field=kTRUE)
Set constant magnetic field and rebuild tracks.
void SetRnrFV(Bool_t x)
Set first-vertex rendering and rebuild tracks.
void LoopToBounds(REveVectorD &p)
Propagate charged particle with momentum p to bounds.
void PrintMagField(Double_t x, Double_t y, Double_t z) const
void SetMaxR(Double_t x)
Set maximum radius and rebuild tracks.
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.
void SetDelta(Double_t x)
Set maximum error and rebuild tracks.
void SetMinAng(Double_t x)
Set maximum step angle and rebuild tracks.
virtual void GoToBounds(REveVectorD &p)
Propagate particle to bounds.
Bool_t LoopToVertex(REveVectorD &v, REveVectorD &p)
Propagate charged particle with momentum p to vertex v.
void SetMagField(Double_t bX, Double_t bY, Double_t bZ)
Set constant magnetic field and rebuild tracks.
Bool_t LineToVertex(REveVectorD &v)
Propagate neutral particle to vertex v.
void SetFitCluster2Ds(Bool_t x)
Set 2D-cluster fitting and rebuild tracks.
REveVector4T A four-vector template without TObject inheritance and virtual functions.
Definition: REveVector.hxx:239
REveVectorT & Sub(const REveVectorT &a, const REveVectorT &b)
Definition: REveVector.hxx:182
TT Dot(const REveVectorT &a) const
Definition: REveVector.hxx:164
Manages Markers.
Definition: TMarker.h:22
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
tbb::task_arena is an alias of tbb::interface7::task_arena, which doesn't allow to forward declare tb...
static constexpr double s
static constexpr double ps
Double_t Sqrt(Double_t x)
Definition: TMath.h:641
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
void UpdateHelix(const REveVectorD &p, const REveVectorD &b, Bool_t full_update, Bool_t enforce_max_step)
Update helix parameters.
void Step(const REveVector4D &v, const REveVectorD &p, REveVector4D &vOut, REveVectorD &pOut)
Step helix for given momentum p from vertex v.
void UpdateRK(const REveVectorD &p, const REveVectorD &b)
Update helix for stepper RungeKutta.
void UpdateCommon(const REveVectorD &p, const REveVectorD &b)
Common update code for helix and RK propagation.