Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TEveTrackPropagator.h
Go to the documentation of this file.
1 // @(#)root/eve:$Id$
2 // Authors: Matevz Tadel & Alja Mrak-Tadel: 2006, 2007
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2007, 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 ROOT_TEveTrackPropagator
13 #define ROOT_TEveTrackPropagator
14 
15 #include "TEveVector.h"
16 #include "TEvePathMark.h"
17 #include "TEveUtil.h"
18 #include "TEveElement.h"
19 #include "TMarker.h"
20 
21 #include <vector>
22 
23 class TEvePointSet;
24 
25 
26 //==============================================================================
27 // TEveMagField
28 //==============================================================================
29 
30 class TEveMagField
31 {
32 protected:
33  Bool_t fFieldConstant;
34 
35 public:
36  TEveMagField() : fFieldConstant(kFALSE) {}
37  virtual ~TEveMagField() {}
38 
39  virtual Bool_t IsConst() const { return fFieldConstant; }
40 
41  virtual void PrintField(Double_t x, Double_t y, Double_t z) const
42  {
43  TEveVector b = GetField(x, y, z);
44  printf("v(%f, %f, %f) B(%f, %f, %f) \n", x, y, z, b.fX, b.fY, b.fZ);
45  }
46 
47  TEveVectorD GetFieldD(const TEveVectorD &v) const { return GetFieldD(v.fX, v.fY, v.fZ); }
48 
49  // Track propgator uses only GetFieldD() and GetMaxFieldMagD(). Have to keep/reuse
50  // GetField() and GetMaxFieldMag() because of backward compatibility.
51 
52  virtual TEveVectorD GetFieldD(Double_t x, Double_t y, Double_t z) const { return GetField(x, y, z); }
53  virtual Double_t GetMaxFieldMagD() const { return GetMaxFieldMag(); } // not abstract because of backward compatibility
54 
55  virtual TEveVector GetField(Float_t, Float_t, Float_t) const { return TEveVector(); }
56  virtual Float_t GetMaxFieldMag() const { return 4; } // not abstract because of backward compatibility
57 
58  ClassDef(TEveMagField, 0); // Abstract interface to magnetic field
59 };
60 
61 
62 //==============================================================================
63 // TEveMagFieldConst
64 //==============================================================================
65 
66 class TEveMagFieldConst : public TEveMagField
67 {
68 protected:
69  TEveVectorD fB;
70 
71 public:
72  TEveMagFieldConst(Double_t x, Double_t y, Double_t z) :
73  TEveMagField(), fB(x, y, z)
74  { fFieldConstant = kTRUE; }
75  virtual ~TEveMagFieldConst() {}
76 
77  using TEveMagField::GetField;
78  virtual TEveVectorD GetFieldD(Double_t /*x*/, Double_t /*y*/, Double_t /*z*/) const { return fB; }
79  virtual Double_t GetMaxFieldMagD() const { return fB.Mag(); };
80 
81  ClassDef(TEveMagFieldConst, 0); // Interface to constant magnetic field.
82 };
83 
84 
85 //==============================================================================
86 // TEveMagFieldDuo
87 //==============================================================================
88 
89 class TEveMagFieldDuo : public TEveMagField
90 {
91 protected:
92  TEveVectorD fBIn;
93  TEveVectorD fBOut;
94  Double_t fR2;
95 
96 public:
97  TEveMagFieldDuo(Double_t r, Double_t bIn, Double_t bOut) :
98  TEveMagField(),
99  fBIn(0,0,bIn), fBOut(0,0,bOut), fR2(r*r)
100  {
101  fFieldConstant = kFALSE;
102  }
103  virtual ~TEveMagFieldDuo() {}
104 
105  using TEveMagField::GetField;
106  virtual TEveVectorD GetFieldD(Double_t x, Double_t y, Double_t /*z*/) const
107  { return ((x*x+y*y)<fR2) ? fBIn : fBOut; }
108 
109  virtual Double_t GetMaxFieldMagD() const
110  { Double_t b1 = fBIn.Mag(), b2 = fBOut.Mag(); return b1 > b2 ? b1 : b2; }
111 
112  ClassDef(TEveMagFieldDuo, 0); // Interface to magnetic field with two different values depending on radius.
113 };
114 
115 
116 //==============================================================================
117 // TEveTrackPropagator
118 //==============================================================================
119 
120 class TEveTrackPropagator : public TEveElementList,
121  public TEveRefBackPtr
122 {
123  friend class TEveTrackPropagatorSubEditor;
124 
125 public:
126  enum EStepper_e { kHelix, kRungeKutta };
127 
128  enum EProjTrackBreaking_e { kPTB_Break, kPTB_UseFirstPointPos, kPTB_UseLastPointPos };
129 
130 protected:
131  struct Helix_t
132  {
133  Int_t fCharge; // Charge of tracked particle.
134  Double_t fMaxAng; // Maximum step angle.
135  Double_t fMaxStep; // Maximum allowed step size.
136  Double_t fDelta; // Maximum error in the middle of the step.
137 
138  Double_t fPhi; // Accumulated angle to check fMaxOrbs by propagator.
139  Bool_t fValid; // Corner case pT~0 or B~0, possible in variable mag field.
140 
141  // ----------------------------------------------------------------
142 
143  // helix parameters
144  Double_t fLam; // Momentum ratio pT/pZ.
145  Double_t fR; // Helix radius in cm.
146  Double_t fPhiStep; // Caluclated from fMinAng and fDelta.
147  Double_t fSin, fCos; // Current sin/cos(phistep).
148 
149  // Runge-Kutta parameters
150  Double_t fRKStep; // Step for Runge-Kutta.
151 
152  // cached
153  TEveVectorD fB; // Current magnetic field, cached.
154  TEveVectorD fE1, fE2, fE3; // Base vectors: E1 -> B dir, E2->pT dir, E3 = E1xE2.
155  TEveVectorD fPt, fPl; // Transverse and longitudinal momentum.
156  Double_t fPtMag; // Magnitude of pT.
157  Double_t fPlMag; // Momentum parallel to mag field.
158  Double_t fLStep; // Transverse step arc-length in cm.
159 
160  // ----------------------------------------------------------------
161 
162  Helix_t();
163 
164  void UpdateCommon(const TEveVectorD & p, const TEveVectorD& b);
165  void UpdateHelix (const TEveVectorD & p, const TEveVectorD& b, Bool_t full_update, Bool_t enforce_max_step);
166  void UpdateRK (const TEveVectorD & p, const TEveVectorD& b);
167 
168  void Step(const TEveVector4D& v, const TEveVectorD& p, TEveVector4D& vOut, TEveVectorD& pOut);
169 
170  Double_t GetStep() { return fLStep * TMath::Sqrt(1 + fLam*fLam); }
171  Double_t GetStep2() { return fLStep * fLStep * (1 + fLam*fLam); }
172  };
173 
174 
175 private:
176  TEveTrackPropagator(const TEveTrackPropagator&); // Not implemented
177  TEveTrackPropagator& operator=(const TEveTrackPropagator&); // Not implemented
178 
179  void DistributeOffset(const TEveVectorD& off, Int_t first_point, Int_t np, TEveVectorD& p);
180 
181 protected:
182  EStepper_e fStepper;
183 
184  TEveMagField* fMagFieldObj;
185  Bool_t fOwnMagFiledObj;
186 
187  // Track extrapolation limits
188  Double_t fMaxR; // Max radius for track extrapolation
189  Double_t fMaxZ; // Max z-coordinate for track extrapolation.
190  Int_t fNMax; // Max steps
191  // Helix limits
192  Double_t fMaxOrbs; // Maximal angular path of tracks' orbits (1 ~ 2Pi).
193 
194  // Path-mark / first-vertex control
195  Bool_t fEditPathMarks; // Show widgets for path-mark control in GUI editor.
196  Bool_t fFitDaughters; // Pass through daughter creation points when extrapolating a track.
197  Bool_t fFitReferences; // Pass through given track-references when extrapolating a track.
198  Bool_t fFitDecay; // Pass through decay point when extrapolating a track.
199  Bool_t fFitCluster2Ds; // Pass through 2D-clusters when extrapolating a track.
200  Bool_t fFitLineSegments; // Pass through line when extrapolating a track.
201  Bool_t fRnrDaughters; // Render daughter path-marks.
202  Bool_t fRnrReferences; // Render track-reference path-marks.
203  Bool_t fRnrDecay; // Render decay path-marks.
204  Bool_t fRnrCluster2Ds; // Render 2D-clusters.
205  Bool_t fRnrFV; // Render first vertex.
206  TMarker fPMAtt; // Marker attributes for rendering of path-marks.
207  TMarker fFVAtt; // Marker attributes for fits vertex.
208 
209  // Handling of discontinuities in projections
210  UChar_t fProjTrackBreaking; // Handling of projected-track breaking.
211  Bool_t fRnrPTBMarkers; // Render break-points on tracks.
212  TMarker fPTBAtt; // Marker attributes for track break-points.
213 
214  // ----------------------------------------------------------------
215 
216  // Propagation, state of current track
217  std::vector<TEveVector4D> fPoints; // Calculated point.
218  std::vector<TEveVector4D> fLastPoints; // Copy of the latest calculated points.
219  TEveVectorD fV; // Start vertex.
220  Helix_t fH; // Helix.
221 
222  void RebuildTracks();
223  void Update(const TEveVector4D& v, const TEveVectorD& p, Bool_t full_update=kFALSE, Bool_t enforce_max_step=kFALSE);
224  void Step(const TEveVector4D &v, const TEveVectorD &p, TEveVector4D &vOut, TEveVectorD &pOut);
225 
226  Bool_t LoopToVertex(TEveVectorD& v, TEveVectorD& p);
227  Bool_t LoopToLineSegment(const TEveVectorD& s, const TEveVectorD& r, TEveVectorD& p);
228  void LoopToBounds(TEveVectorD& p);
229 
230  Bool_t LineToVertex (TEveVectorD& v);
231  void LineToBounds (TEveVectorD& p);
232 
233  void StepRungeKutta(Double_t step, Double_t* vect, Double_t* vout);
234 
235  Bool_t HelixIntersectPlane(const TEveVectorD& p, const TEveVectorD& point, const TEveVectorD& normal,
236  TEveVectorD&itsect);
237  Bool_t LineIntersectPlane(const TEveVectorD& p, const TEveVectorD& point, const TEveVectorD& normal,
238  TEveVectorD& itsect);
239  Bool_t PointOverVertex(const TEveVector4D& v0, const TEveVector4D& v, Double_t* p=0);
240 
241  void ClosestPointFromVertexToLineSegment(const TEveVectorD& v, const TEveVectorD& s, const TEveVectorD& r, Double_t rMagInv, TEveVectorD& c);
242  Bool_t ClosestPointBetweenLines(const TEveVectorD&, const TEveVectorD&, const TEveVectorD&, const TEveVectorD&, TEveVectorD& out);
243 
244 public:
245  TEveTrackPropagator(const char* n="TEveTrackPropagator", const char* t="",
246  TEveMagField* field=0, Bool_t own_field=kTRUE);
247  virtual ~TEveTrackPropagator();
248 
249  virtual void OnZeroRefCount();
250 
251  virtual void CheckReferenceCount(const TEveException& eh="TEveElement::CheckReferenceCount ");
252 
253  virtual void ElementChanged(Bool_t update_scenes=kTRUE, Bool_t redraw=kFALSE);
254 
255  // propagation
256  void InitTrack(const TEveVectorD& v, Int_t charge);
257  void ResetTrack();
258 
259  Int_t GetCurrentPoint() const;
260  Double_t GetTrackLength(Int_t start_point=0, Int_t end_point=-1) const;
261 
262  virtual void GoToBounds(TEveVectorD& p);
263  virtual Bool_t GoToVertex(TEveVectorD& v, TEveVectorD& p);
264  virtual Bool_t GoToLineSegment(const TEveVectorD& s, const TEveVectorD& r, TEveVectorD& p);
265 
266  // TEveVectorF wrappers
267  void InitTrack(const TEveVectorF& v, Int_t charge);
268  void GoToBounds(TEveVectorF& p);
269  Bool_t GoToVertex(TEveVectorF& v, TEveVectorF&p);
270  Bool_t GoToLineSegment(const TEveVectorF& s, const TEveVectorF& r, TEveVectorF& p);
271 
272  Bool_t IntersectPlane(const TEveVectorD& p, const TEveVectorD& point, const TEveVectorD& normal,
273  TEveVectorD& itsect);
274 
275  void FillPointSet(TEvePointSet* ps) const;
276 
277  void SetStepper(EStepper_e s) { fStepper = s; }
278 
279  void SetMagField(Double_t bX, Double_t bY, Double_t bZ);
280  void SetMagField(Double_t b) { SetMagField(0, 0, b); }
281  void SetMagFieldObj(TEveMagField* field, Bool_t own_field=kTRUE);
282 
283  void SetMaxR(Double_t x);
284  void SetMaxZ(Double_t x);
285  void SetMaxOrbs(Double_t x);
286  void SetMinAng(Double_t x);
287  void SetMaxAng(Double_t x);
288  void SetMaxStep(Double_t x);
289  void SetDelta(Double_t x);
290 
291  void SetEditPathMarks(Bool_t x) { fEditPathMarks = x; }
292  void SetRnrDaughters(Bool_t x);
293  void SetRnrReferences(Bool_t x);
294  void SetRnrDecay(Bool_t x);
295  void SetRnrCluster2Ds(Bool_t x);
296  void SetFitDaughters(Bool_t x);
297  void SetFitReferences(Bool_t x);
298  void SetFitDecay(Bool_t x);
299  void SetFitCluster2Ds(Bool_t x);
300  void SetFitLineSegments(Bool_t x);
301  void SetRnrFV(Bool_t x);
302  void SetProjTrackBreaking(UChar_t x);
303  void SetRnrPTBMarkers(Bool_t x);
304 
305  TEveVectorD GetMagField(Double_t x, Double_t y, Double_t z) { return fMagFieldObj->GetField(x, y, z); }
306  void PrintMagField(Double_t x, Double_t y, Double_t z) const;
307 
308  EStepper_e GetStepper() const { return fStepper;}
309 
310  Double_t GetMaxR() const { return fMaxR; }
311  Double_t GetMaxZ() const { return fMaxZ; }
312  Double_t GetMaxOrbs() const { return fMaxOrbs; }
313  Double_t GetMinAng() const;
314  Double_t GetMaxAng() const { return fH.fMaxAng; }
315  Double_t GetMaxStep() const { return fH.fMaxStep; }
316  Double_t GetDelta() const { return fH.fDelta; }
317 
318  Bool_t GetEditPathMarks() const { return fEditPathMarks; }
319  Bool_t GetRnrDaughters() const { return fRnrDaughters; }
320  Bool_t GetRnrReferences() const { return fRnrReferences; }
321  Bool_t GetRnrDecay() const { return fRnrDecay; }
322  Bool_t GetRnrCluster2Ds() const { return fRnrCluster2Ds; }
323  Bool_t GetFitDaughters() const { return fFitDaughters; }
324  Bool_t GetFitReferences() const { return fFitReferences; }
325  Bool_t GetFitDecay() const { return fFitDecay; }
326  Bool_t GetFitCluster2Ds() const { return fFitCluster2Ds; }
327  Bool_t GetFitLineSegments() const { return fFitLineSegments; }
328  Bool_t GetRnrFV() const { return fRnrFV; }
329  UChar_t GetProjTrackBreaking() const { return fProjTrackBreaking; }
330  Bool_t GetRnrPTBMarkers() const { return fRnrPTBMarkers; }
331 
332  TMarker& RefPMAtt() { return fPMAtt; }
333  TMarker& RefFVAtt() { return fFVAtt; }
334  TMarker& RefPTBAtt() { return fPTBAtt; }
335 
336  const std::vector<TEveVector4D>& GetLastPoints() const { return fLastPoints; }
337 
338  static Bool_t IsOutsideBounds(const TEveVectorD& point, Double_t maxRsqr, Double_t maxZ);
339 
340  static Double_t fgDefMagField; // Default value for constant solenoid magnetic field.
341  static const Double_t fgkB2C; // Constant for conversion of momentum to curvature.
342  static TEveTrackPropagator fgDefault; // Default track propagator.
343 
344  static Double_t fgEditorMaxR; // Max R that can be set in GUI editor.
345  static Double_t fgEditorMaxZ; // Max Z that can be set in GUI editor.
346 
347  ClassDef(TEveTrackPropagator, 0); // Calculates path of a particle taking into account special path-marks and imposed boundaries.
348 };
349 
350 //______________________________________________________________________________
351 inline Bool_t TEveTrackPropagator::IsOutsideBounds(const TEveVectorD& point,
352  Double_t maxRsqr,
353  Double_t maxZ)
354 {
355  // Return true if point% is outside of cylindrical bounds detrmined by
356  // square radius and z.
357 
358  return TMath::Abs(point.fZ) > maxZ ||
359  point.fX*point.fX + point.fY*point.fY > maxRsqr;
360 }
361 
362 //______________________________________________________________________________
363 inline Bool_t TEveTrackPropagator::PointOverVertex(const TEveVector4D &v0,
364  const TEveVector4D &v,
365  Double_t *p)
366 {
367  static const Double_t kMinPl = 1e-5;
368 
369  TEveVectorD dv; dv.Sub(v0, v);
370 
371  Double_t dotV;
372 
373  if (TMath::Abs(fH.fPlMag) > kMinPl)
374  {
375  // Use longitudinal momentum to determine crossing point.
376  // Works ok for spiraling helices, also for loopers.
377 
378  dotV = fH.fE1.Dot(dv);
379  if (fH.fPlMag < 0)
380  dotV = -dotV;
381  }
382  else
383  {
384  // Use full momentum, which is pT, under this conditions.
385 
386  dotV = fH.fE2.Dot(dv);
387  }
388 
389  if (p)
390  *p = dotV;
391 
392  return dotV < 0;
393 }
394 
395 #endif