Logo ROOT  
Reference Guide
TEveTrackPropagator.cxx
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
13#include "TEveTrackPropagator.h"
14#include "TEveTrack.h"
15#include "TEveTrans.h"
16
17#include "TMath.h"
18
19#include <cassert>
20
21/** \class TEveMagField
22\ingroup TEve
23Abstract base-class for interfacing to magnetic field needed by the
24TEveTrackPropagator.
25
26To implement your own version, redefine the following virtual functions:
27 virtual Double_t GetMaxFieldMagD() const;
28 virtual TEveVectorD GetFieldD(Double_t x, Double_t y, Double_t z) const;
29
30See sub-classes TEveMagFieldConst and TEveMagFieldDuo for two simple implementations.
31
32NOTE: For backward compatibility float versions are the primary
33sources of field information in this base-class. The float versions are
34not used in TEve and can be ignored in sub-classes.
35
36NOTE: Magnetic field direction convention is inverted.
37*/
38
40
41/** \class TEveMagFieldConst
42\ingroup TEve
43Implements constant magnetic field, given by a vector fB.
44
45NOTE: Magnetic field direction convention is inverted.
46*/
47
49
50/** \class TEveMagFieldDuo
51\ingroup TEve
52Implements constant magnetic filed that switches on given axial radius fR2
53from vector fBIn to fBOut.
54
55NOTE: Magnetic field direction convention is inverted.
56*/
57
59
60namespace
61{
62 //const Double_t kBMin = 1e-6;
63 const Double_t kPtMinSqr = 1e-20;
64 const Double_t kAMin = 1e-10;
65 const Double_t kStepEps = 1e-3;
66}
67
68////////////////////////////////////////////////////////////////////////////////
69/// Default constructor.
70
72 fCharge(0),
73 fMaxAng(45), fMaxStep(20.f), fDelta(0.1),
74 fPhi(0), fValid(kFALSE),
75 fLam(-1), fR(-1), fPhiStep(-1), fSin(-1), fCos(-1),
76 fRKStep(20.0),
77 fPtMag(-1), fPlMag(-1), fLStep(-1)
78{
79}
80
81////////////////////////////////////////////////////////////////////////////////
82/// Common update code for helix and RK propagation.
83
85{
86 fB = b;
87
88 // base vectors
89 fE1 = b;
90 fE1.Normalize();
91 fPlMag = p.Dot(fE1);
92 fPl = fE1*fPlMag;
93
94 fPt = p - fPl;
95 fPtMag = fPt.Mag();
96 fE2 = fPt;
97 fE2.Normalize();
98}
99
100////////////////////////////////////////////////////////////////////////////////
101/// Update helix parameters.
102
104 Bool_t full_update, Bool_t enforce_max_step)
105{
106 UpdateCommon(p, b);
107
108 // helix parameters
109 TMath::Cross(fE1.Arr(), fE2.Arr(), fE3.Arr());
110 if (fCharge < 0) fE3.NegateXYZ();
111
112 if (full_update)
113 {
114 using namespace TMath;
115 Double_t a = fgkB2C * b.Mag() * Abs(fCharge);
116 if (a > kAMin && fPtMag*fPtMag > kPtMinSqr)
117 {
118 fValid = kTRUE;
119
120 fR = Abs(fPtMag / a);
121 fLam = fPlMag / fPtMag;
122
123 // get phi step, compare fMaxAng with fDelta
124 fPhiStep = fMaxAng * DegToRad();
125 if (fR > fDelta)
126 {
127 Double_t ang = 2.0 * ACos(1.0f - fDelta/fR);
128 if (ang < fPhiStep)
129 fPhiStep = ang;
130 }
131
132 // check max step size
133 Double_t curr_step = fR * fPhiStep * Sqrt(1.0f + fLam*fLam);
134 if (curr_step > fMaxStep || enforce_max_step)
135 fPhiStep *= fMaxStep / curr_step;
136
137 fLStep = fR * fPhiStep * fLam;
138 fSin = Sin(fPhiStep);
139 fCos = Cos(fPhiStep);
140 }
141 else
142 {
143 fValid = kFALSE;
144 }
145 }
146}
147
148////////////////////////////////////////////////////////////////////////////////
149/// Update helix for stepper RungeKutta.
150
152{
153 UpdateCommon(p, b);
154
155 if (fCharge)
156 {
157 fValid = kTRUE;
158
159 // cached values for propagator
160 fB = b;
161 fPlMag = p.Dot(fB);
162 }
163 else
164 {
165 fValid = kFALSE;
166 }
167}
168
169////////////////////////////////////////////////////////////////////////////////
170/// Step helix for given momentum p from vertex v.
171
173 TEveVector4D& vOut, TEveVectorD& pOut)
174{
175 vOut = v;
176
177 if (fValid)
178 {
179 TEveVectorD d = fE2*(fR*fSin) + fE3*(fR*(1-fCos)) + fE1*fLStep;
180 vOut += d;
181 vOut.fT += TMath::Abs(fLStep);
182
183 pOut = fPl + fE2*(fPtMag*fCos) + fE3*(fPtMag*fSin);
184
185 fPhi += fPhiStep;
186 }
187 else
188 {
189 // case: pT < kPtMinSqr or B < kBMin
190 // might happen if field directon changes pT ~ 0 or B becomes zero
191 vOut += p * (fMaxStep / p.Mag());
192 vOut.fT += fMaxStep;
193 pOut = p;
194 }
195}
196
197/** \class TEveTrackPropagator
198\ingroup TEve
199Holding structure for a number of track rendering parameters.
200Calculates path taking into account the parameters.
201
202NOTE: Magnetic field direction convention is inverted.
203
204This is decoupled from TEveTrack/TEveTrackList to allow sharing of the
205Propagator among several instances. Back references are kept so the tracks
206can be recreated when the parameters change.
207
208TEveTrackList has Get/Set methods for RnrStlye. TEveTrackEditor and
209TEveTrackListEditor provide editor access.
210
211Enum EProjTrackBreaking_e and member fProjTrackBreaking specify whether 2D
212projected tracks get broken into several segments when the projected space
213consists of separate domains (like Rho-Z). The track-breaking is enabled by
214default.
215*/
216
218
220const Double_t TEveTrackPropagator::fgkB2C = 0.299792458e-2;
222
225
226////////////////////////////////////////////////////////////////////////////////
227/// Default constructor.
228
230 TEveMagField *field, Bool_t own_field) :
231 TEveElementList(n, t),
233
235 fMagFieldObj(field),
236 fOwnMagFiledObj(own_field),
237
238 fMaxR (350), fMaxZ (450),
239 fNMax (4096), fMaxOrbs (0.5),
240
247 fRnrFV (kFALSE),
248 fPMAtt(), fFVAtt(),
249
251
252 fV()
253{
257
261
265
266 if (fMagFieldObj == 0) {
269 }
270}
271
272////////////////////////////////////////////////////////////////////////////////
273/// Destructor.
274
276{
277 if (fOwnMagFiledObj)
278 {
279 delete fMagFieldObj;
280 }
281}
282
283////////////////////////////////////////////////////////////////////////////////
284/// Virtual from TEveRefBackPtr - track reference count has reached zero.
285
287{
288 CheckReferenceCount("TEveTrackPropagator::OnZeroRefCount ");
289}
290
291////////////////////////////////////////////////////////////////////////////////
292/// Check reference count - virtual from TEveElement.
293/// Must also take into account references from TEveRefBackPtr.
294
296{
297 if (fRefCount <= 0)
298 {
300 }
301}
302
303////////////////////////////////////////////////////////////////////////////////
304/// Element-change notification.
305/// Stamp all tracks as requiring display-list regeneration.
306/// Virtual from TEveElement.
307
309{
310 TEveTrack* track;
311 RefMap_i i = fBackRefs.begin();
312 while (i != fBackRefs.end())
313 {
314 track = dynamic_cast<TEveTrack*>(i->first);
315 track->StampObjProps();
316 ++i;
317 }
318 TEveElementList::ElementChanged(update_scenes, redraw);
319}
320
321////////////////////////////////////////////////////////////////////////////////
322/// Initialize internal data-members for given particle parameters.
323
325{
326 fV = v;
327 fPoints.push_back(fV);
328
329 // init helix
330 fH.fPhi = 0;
331 fH.fCharge = charge;
332}
333
334////////////////////////////////////////////////////////////////////////////////
335/// TEveVectorF wrapper.
336
338{
339 TEveVectorD vd(v);
340 InitTrack(vd, charge);
341}
342
343////////////////////////////////////////////////////////////////////////////////
344/// Reset cache holding particle trajectory.
345
347{
348 fLastPoints.clear();
349 fPoints.swap(fLastPoints);
350
351 // reset helix
352 fH.fPhi = 0;
353}
354
355////////////////////////////////////////////////////////////////////////////////
356/// Get index of current point on track.
357
359{
360 return fPoints.size() - 1;
361}
362
363////////////////////////////////////////////////////////////////////////////////
364/// Calculate track length from start_point to end_point.
365/// If end_point is less than 0, distance to the end is returned.
366
368{
369 if (end_point < 0) end_point = fPoints.size() - 1;
370
371 Double_t sum = 0;
372 for (Int_t i = start_point; i < end_point; ++i)
373 {
374 sum += (fPoints[i+1] - fPoints[i]).Mag();
375 }
376 return sum;
377}
378
379////////////////////////////////////////////////////////////////////////////////
380/// Propagate particle with momentum p to vertex v.
381
383{
384 Update(fV, p, kTRUE);
385
386 if ((v-fV).Mag() < kStepEps)
387 {
388 fPoints.push_back(v);
389 return kTRUE;
390 }
391
392 return fH.fValid ? LoopToVertex(v, p) : LineToVertex(v);
393}
394
395////////////////////////////////////////////////////////////////////////////////
396/// Propagate particle with momentum p to line with start point s and vector r
397/// to the second point.
398
400{
401 Update(fV, p, kTRUE);
402
403 if (!fH.fValid)
404 {
408 return kTRUE;
409 }
410 else
411 {
412 return LoopToLineSegment(s, r, p);
413 }
414}
415
416////////////////////////////////////////////////////////////////////////////////
417/// TEveVectorF wrapper.
418
420{
421 TEveVectorD vd(v), pd(p);
422 Bool_t result = GoToVertex(vd, pd);
423 v = vd; p = pd;
424 return result;
425}
426
427////////////////////////////////////////////////////////////////////////////////
428/// TEveVectorF wrapper.
429
431{
432 TEveVectorD sd(s), rd(r), pd(p);
433 Bool_t result = GoToLineSegment(sd, rd, pd);
434 p = pd;
435 return result;
436}
437
438////////////////////////////////////////////////////////////////////////////////
439/// Propagate particle to bounds.
440/// Return TRUE if hit bounds.
441
443{
444 Update(fV, p, kTRUE);
445
447}
448
449////////////////////////////////////////////////////////////////////////////////
450/// TEveVectorF wrapper.
451
453{
454 TEveVectorD pd(p);
455 GoToBounds(pd);
456 p = pd;
457}
458
459////////////////////////////////////////////////////////////////////////////////
460/// Update helix / B-field projection state.
461
463 Bool_t full_update, Bool_t enforce_max_step)
464{
465 if (fStepper == kHelix)
466 {
467 fH.UpdateHelix(p, fMagFieldObj->GetFieldD(v), !fMagFieldObj->IsConst() || full_update, enforce_max_step);
468 }
469 else
470 {
472
473 if (full_update)
474 {
475 using namespace TMath;
476
478 if (a > kAMin)
479 {
480 fH.fR = p.Mag() / a;
481
482 // get phi step, compare fDelta with MaxAng
484 if (fH.fR > fH.fDelta )
485 {
486 Double_t ang = 2.0 * ACos(1.0f - fH.fDelta/fH.fR);
487 if (ang < fH.fPhiStep)
488 fH.fPhiStep = ang;
489 }
490
491 // check against maximum step-size
492 fH.fRKStep = fH.fR * fH.fPhiStep * Sqrt(1 + fH.fLam*fH.fLam);
493 if (fH.fRKStep > fH.fMaxStep || enforce_max_step)
494 {
497 }
498 }
499 else
500 {
502 }
503 }
504 }
505}
506
507////////////////////////////////////////////////////////////////////////////////
508/// Wrapper to step helix.
509
511{
512 if (fStepper == kHelix)
513 {
514 fH.Step(v, p, vOut, pOut);
515 }
516 else
517 {
518 Double_t vecRKIn[7];
519 vecRKIn[0] = v.fX;
520 vecRKIn[1] = v.fY;
521 vecRKIn[2] = v.fZ;
522 Double_t pm = p.Mag();
523 Double_t nm = 1.0 / pm;
524 vecRKIn[3] = p.fX*nm;
525 vecRKIn[4] = p.fY*nm;
526 vecRKIn[5] = p.fZ*nm;
527 vecRKIn[6] = p.Mag();
528
529 Double_t vecRKOut[7];
530 StepRungeKutta(fH.fRKStep, vecRKIn, vecRKOut);
531
532 vOut.fX = vecRKOut[0];
533 vOut.fY = vecRKOut[1];
534 vOut.fZ = vecRKOut[2];
535 vOut.fT = v.fT + fH.fRKStep;
536 pm = vecRKOut[6];
537 pOut.fX = vecRKOut[3]*pm;
538 pOut.fY = vecRKOut[4]*pm;
539 pOut.fZ = vecRKOut[5]*pm;
540 }
541}
542
543////////////////////////////////////////////////////////////////////////////////
544/// Propagate charged particle with momentum p to bounds.
545/// It is expected that Update() with full-update was called before.
546
548{
549 const Double_t maxRsq = fMaxR*fMaxR;
550
551 TEveVector4D currV(fV);
552 TEveVector4D forwV(fV);
553 TEveVectorD forwP (p);
554
555 Int_t np = fPoints.size();
556 Double_t maxPhi = fMaxOrbs*TMath::TwoPi();
557
558 while (fH.fPhi < maxPhi && np<fNMax)
559 {
560 Step(currV, p, forwV, forwP);
561
562 // cross R
563 if (forwV.Perp2() > maxRsq)
564 {
565 Float_t t = (fMaxR - currV.R()) / (forwV.R() - currV.R());
566 if (t < 0 || t > 1)
567 {
568 Warning("HelixToBounds", "In MaxR crossing expected t>=0 && t<=1: t=%f, r1=%f, r2=%f, MaxR=%f.",
569 t, currV.R(), forwV.R(), fMaxR);
570 return;
571 }
572 TEveVectorD d(forwV);
573 d -= currV;
574 d *= t;
575 d += currV;
576 fPoints.push_back(d);
577 return;
578 }
579
580 // cross Z
581 else if (TMath::Abs(forwV.fZ) > fMaxZ)
582 {
583 Double_t t = (fMaxZ - TMath::Abs(currV.fZ)) / TMath::Abs((forwV.fZ - currV.fZ));
584 if (t < 0 || t > 1)
585 {
586 Warning("HelixToBounds", "In MaxZ crossing expected t>=0 && t<=1: t=%f, z1=%f, z2=%f, MaxZ=%f.",
587 t, currV.fZ, forwV.fZ, fMaxZ);
588 return;
589 }
590 TEveVectorD d(forwV -currV);
591 d *= t;
592 d += currV;
593 fPoints.push_back(d);
594 return;
595 }
596
597 currV = forwV;
598 p = forwP;
599 Update(currV, p);
600
601 fPoints.push_back(currV);
602 ++np;
603 }
604}
605
606////////////////////////////////////////////////////////////////////////////////
607/// Propagate charged particle with momentum p to vertex v.
608/// It is expected that Update() with full-update was called before.
609
611{
612 const Double_t maxRsq = fMaxR * fMaxR;
613
614 TEveVector4D currV(fV);
615 TEveVector4D forwV(fV);
616 TEveVectorD forwP(p);
617
618 Int_t first_point = fPoints.size();
619 Int_t np = first_point;
620
621 Double_t prod0=0, prod1;
622
623 do
624 {
625 Step(currV, p, forwV, forwP);
626 Update(forwV, forwP);
627
628 if (PointOverVertex(v, forwV, &prod1))
629 {
630 break;
631 }
632
633 if (IsOutsideBounds(forwV, maxRsq, fMaxZ))
634 {
635 fV = currV;
636 return kFALSE;
637 }
638
639 fPoints.push_back(forwV);
640 currV = forwV;
641 p = forwP;
642 prod0 = prod1;
643 ++np;
644 } while (np < fNMax);
645
646 // make the remaining fractional step
647 if (np > first_point)
648 {
649 if ((v - currV).Mag() > kStepEps)
650 {
651 Double_t step_frac = prod0 / (prod0 - prod1);
652 if (step_frac > 0)
653 {
654 // Step for fraction of previous step size.
655 // We pass 'enforce_max_step' flag to Update().
656 Float_t orig_max_step = fH.fMaxStep;
657 fH.fMaxStep = step_frac * (forwV - currV).Mag();
658 Update(currV, p, kTRUE, kTRUE);
659 Step(currV, p, forwV, forwP);
660 p = forwP;
661 currV = forwV;
662 fPoints.push_back(currV);
663 ++np;
664 fH.fMaxStep = orig_max_step;
665 }
666
667 // Distribute offset to desired crossing point over all segment.
668
669 TEveVectorD off(v - currV);
670 off *= 1.0f / currV.fT;
671 DistributeOffset(off, first_point, np, p);
672 fV = v;
673 return kTRUE;
674 }
675 }
676
677 fPoints.push_back(v);
678 fV = v;
679 return kTRUE;
680}
681
682////////////////////////////////////////////////////////////////////////////////
683/// Propagate charged particle with momentum p to line segment with point s and
684/// vector r to the second point. It is expected that Update() with full-update
685/// was called before. Returns kFALSE if hits bounds.
686
688{
689 const Double_t maxRsq = fMaxR * fMaxR;
690 const Double_t rMagInv = 1./r.Mag();
691
692 TEveVector4D currV(fV);
693 TEveVector4D forwV(fV);
694 TEveVectorD forwP(p);
695
696 Int_t first_point = fPoints.size();
697 Int_t np = first_point;
698
699 TEveVectorD forwC;
700 TEveVectorD currC;
701 do
702 {
703 Step(currV, p, forwV, forwP);
704 Update(forwV, forwP);
705
706 ClosestPointFromVertexToLineSegment(forwV, s, r, rMagInv, forwC);
707
708 // check forwV is over segment with orthogonal component of
709 // momentum to vector r
710 TEveVectorD b = r; b.Normalize();
711 Double_t x = forwP.Dot(b);
712 TEveVectorD pTPM = forwP - x*b;
713 if (pTPM.Dot(forwC - forwV) < 0)
714 {
715 break;
716 }
717
718 if (IsOutsideBounds(forwV, maxRsq, fMaxZ))
719 {
720 fV = currV;
721 return kFALSE;
722 }
723
724 fPoints.push_back(forwV);
725 currV = forwV;
726 p = forwP;
727 currC = forwC;
728 ++np;
729 } while (np < fNMax);
730
731 // Get closest point on segment relative to line with forw and currV points.
733 ClosestPointBetweenLines(s, r, currV, forwV - currV, v);
734
735 // make the remaining fractional step
736 if (np > first_point)
737 {
738 if ((v - currV).Mag() > kStepEps)
739 {
740 TEveVector last_step = forwV - currV;
741 TEveVector delta = v - currV;
742 Double_t step_frac = last_step.Dot(delta) / last_step.Mag2();
743 if (step_frac > 0)
744 {
745 // Step for fraction of previous step size.
746 // We pass 'enforce_max_step' flag to Update().
747 Float_t orig_max_step = fH.fMaxStep;
748 fH.fMaxStep = step_frac * (forwV - currV).Mag();
749 Update(currV, p, kTRUE, kTRUE);
750 Step(currV, p, forwV, forwP);
751 p = forwP;
752 currV = forwV;
753 fPoints.push_back(currV);
754 ++np;
755 fH.fMaxStep = orig_max_step;
756 }
757
758 // Distribute offset to desired crossing point over all segment.
759
760 TEveVectorD off(v - currV);
761 off *= 1.0f / currV.fT;
762 DistributeOffset(off, first_point, np, p);
763 fV = v;
764 return kTRUE;
765 }
766 }
767
768 fPoints.push_back(v);
769 fV = v;
770 return kTRUE;
771}
772
773////////////////////////////////////////////////////////////////////////////////
774/// Distribute offset between first and last point index and rotate
775/// momentum.
776
778{
779 // Calculate the required momentum rotation.
780 // lpd - last-points-delta
781 TEveVectorD lpd0(fPoints[np-1]);
782 lpd0 -= fPoints[np-2];
783 lpd0.Normalize();
784
785 for (Int_t i = first_point; i < np; ++i)
786 {
787 fPoints[i] += off * fPoints[i].fT;
788 }
789
790 TEveVectorD lpd1(fPoints[np-1]);
791 lpd1 -= fPoints[np-2];
792 lpd1.Normalize();
793
795 tt.SetupFromToVec(lpd0, lpd1);
796
797 // TEveVectorD pb4(p);
798 // printf("Rotating momentum: p0 = "); p.Dump();
799 tt.RotateIP(p);
800 // printf(" p1 = "); p.Dump();
801 // printf(" n1=%f, n2=%f, dp = %f deg\n", pb4.Mag(), p.Mag(),
802 // TMath::RadToDeg()*TMath::ACos(p.Dot(pb4)/(pb4.Mag()*p.Mag())));
803}
804
805////////////////////////////////////////////////////////////////////////////////
806/// Propagate neutral particle to vertex v.
807
809{
810 TEveVector4D currV = v;
811
812 currV.fX = v.fX;
813 currV.fY = v.fY;
814 currV.fZ = v.fZ;
815 fPoints.push_back(currV);
816
817 fV = v;
818 return kTRUE;
819}
820
821////////////////////////////////////////////////////////////////////////////////
822/// Propagate neutral particle with momentum p to bounds.
823
825{
826 Double_t tZ = 0, tR = 0, tB = 0;
827
828 // time where particle intersect +/- fMaxZ
829 if (p.fZ > 0)
830 tZ = (fMaxZ - fV.fZ) / p.fZ;
831 else if (p.fZ < 0)
832 tZ = - (fMaxZ + fV.fZ) / p.fZ;
833 else
834 tZ = 1e99;
835
836 // time where particle intersects cylinder
837 Double_t a = p.fX*p.fX + p.fY*p.fY;
838 Double_t b = 2.0 * (fV.fX*p.fX + fV.fY*p.fY);
840 Double_t d = b*b - 4.0*a*c;
841 if (d >= 0) {
842 Double_t sqrtD = TMath::Sqrt(d);
843 tR = (-b - sqrtD) / (2.0 * a);
844 if (tR < 0) {
845 tR = (-b + sqrtD) / (2.0 * a);
846 }
847 tB = tR < tZ ? tR : tZ; // compare the two times
848 } else {
849 tB = tZ;
850 }
851 TEveVectorD nv(fV.fX + p.fX*tB, fV.fY + p.fY*tB, fV.fZ + p.fZ*tB);
852 LineToVertex(nv);
853}
854
855////////////////////////////////////////////////////////////////////////////////
856/// Intersect helix with a plane. Current position and argument p define
857/// the helix.
858
860 const TEveVectorD& point,
861 const TEveVectorD& normal,
862 TEveVectorD& itsect)
863{
864 TEveVectorD pos(fV);
865 TEveVectorD mom(p);
866 if (fMagFieldObj->IsConst())
868
869 TEveVectorD n(normal);
870 TEveVectorD delta = pos - point;
871 Double_t d = delta.Dot(n);
872 if (d > 0) {
873 n.NegateXYZ(); // Turn normal around so that we approach from negative side of the plane
874 d = -d;
875 }
876
877 TEveVector4D forwV;
878 TEveVectorD forwP;
879 TEveVector4D pos4(pos);
880 while (kTRUE)
881 {
882 Update(pos4, mom);
883 Step(pos4, mom, forwV , forwP);
884 Double_t new_d = (forwV - point).Dot(n);
885 if (new_d < d)
886 {
887 // We are going further away ... fail intersect.
888 Warning("HelixIntersectPlane", "going away from the plane.");
889 return kFALSE;
890 }
891 if (new_d > 0)
892 {
893 delta = forwV - pos4;
894 itsect = pos4 + delta * ((point - pos4).Dot(n) / delta.Dot(n));
895
896 return kTRUE;
897 }
898 pos4 = forwV;
899 mom = forwP;
900 }
901}
902
903////////////////////////////////////////////////////////////////////////////////
904/// Intersect line with a plane. Current position and argument p define
905/// the line.
906
908 const TEveVectorD& point,
909 const TEveVectorD& normal,
910 TEveVectorD& itsect)
911{
912 TEveVectorD pos(fV.fX, fV.fY, fV.fZ);
913 TEveVectorD delta = point - pos;
914
915 Double_t pn = p.Dot(normal);
916 if (pn == 0)
917 {
918 return kFALSE;
919 }
920 Double_t t = delta.Dot(normal) / pn;
921 if (t < 0) {
922 return kFALSE;
923 } else {
924 itsect = pos + p*t;
925 return kTRUE;
926 }
927}
928
929////////////////////////////////////////////////////////////////////////////////
930/// Find intersection of currently propagated track with a plane.
931/// Current track position is used as starting point.
932///
933/// Args:
934/// - p - track momentum to use for extrapolation
935/// - point - a point on a plane
936/// - normal - normal of the plane
937/// - itsect - output, point of intersection
938/// Returns:
939/// - kFALSE if intersection can not be found, kTRUE otherwise.
940
942 const TEveVectorD& point,
943 const TEveVectorD& normal,
944 TEveVectorD& itsect)
945{
946 if (fH.fCharge && fMagFieldObj)
947 return HelixIntersectPlane(p, point, normal, itsect);
948 else
949 return LineIntersectPlane(p, point, normal, itsect);
950}
951
952////////////////////////////////////////////////////////////////////////////////
953/// Get closest point from given vertex v to line segment defined with s and r.
954/// Argument rMagInv is cached. rMagInv= 1./rMag()
955
957 const TEveVectorD& s,
958 const TEveVectorD& r,
959 Double_t rMagInv,
960 TEveVectorD& c)
961{
962 TEveVectorD dir = v - s;
963 TEveVectorD b1 = r * rMagInv;
964
965 // parallel distance
966 Double_t dot = dir.Dot(b1);
967 TEveVectorD dirI = dot * b1;
968
969 Double_t facX = dot * rMagInv;
970
971 if (facX <= 0)
972 c = s;
973 else if (facX >= 1)
974 c = s + r;
975 else
976 c = s + dirI;
977}
978
979////////////////////////////////////////////////////////////////////////////////
980/// Get closest point on line defined with vector p0 and u.
981/// Return false if the point is forced on the line segment.
982
984 const TEveVectorD& u,
985 const TEveVectorD& q0,
986 const TEveVectorD& v,
987 TEveVectorD& out)
988{
989 TEveVectorD w0 = p0 -q0;
990 Double_t a = u.Mag2();
991 Double_t b = u.Dot(v);
992 Double_t c = v.Mag2();
993 Double_t d = u.Dot(w0);
994 Double_t e = v.Dot(w0);
995
996 Double_t x = (b*e - c*d)/(a*c -b*b);
997 Bool_t force = (x < 0 || x > 1);
998 out = p0 + TMath::Range(0., 1., x) * u;
999 return force;
1000}
1001
1002////////////////////////////////////////////////////////////////////////////////
1003/// Reset ps and populate it with points in propagation cache.
1004
1006{
1007 Int_t size = TMath::Min(fNMax, (Int_t)fPoints.size());
1008 ps->Reset(size);
1009 for (Int_t i = 0; i < size; ++i)
1010 {
1011 const TEveVector4D& v = fPoints[i];
1012 ps->SetNextPoint(v.fX, v.fY, v.fZ);
1013 }
1014}
1015
1016////////////////////////////////////////////////////////////////////////////////
1017/// Rebuild all tracks using this render-style.
1018
1020{
1021 TEveTrack* track;
1022 RefMap_i i = fBackRefs.begin();
1023 while (i != fBackRefs.end())
1024 {
1025 track = dynamic_cast<TEveTrack*>(i->first);
1026 track->MakeTrack();
1027 track->StampObjProps();
1028 ++i;
1029 }
1030}
1031
1032////////////////////////////////////////////////////////////////////////////////
1033/// Set constant magnetic field and rebuild tracks.
1034
1036{
1037 SetMagFieldObj(new TEveMagFieldConst(bX, bY, bZ));
1038}
1039
1040////////////////////////////////////////////////////////////////////////////////
1041/// Set constant magnetic field and rebuild tracks.
1042
1044{
1046
1047 fMagFieldObj = field;
1048 fOwnMagFiledObj = own_field;
1049
1050 RebuildTracks();
1051}
1052
1053////////////////////////////////////////////////////////////////////////////////
1054
1056{
1058}
1059
1060////////////////////////////////////////////////////////////////////////////////
1061/// Set maximum radius and rebuild tracks.
1062
1064{
1065 fMaxR = x;
1066 RebuildTracks();
1067}
1068
1069////////////////////////////////////////////////////////////////////////////////
1070/// Set maximum z and rebuild tracks.
1071
1073{
1074 fMaxZ = x;
1075 RebuildTracks();
1076}
1077
1078////////////////////////////////////////////////////////////////////////////////
1079/// Set maximum number of orbits and rebuild tracks.
1080
1082{
1083 fMaxOrbs = x;
1084 RebuildTracks();
1085}
1086
1087////////////////////////////////////////////////////////////////////////////////
1088/// Set maximum step angle and rebuild tracks.
1089/// WARNING -- this method / variable was mis-named.
1090
1092{
1093 Warning("SetMinAng", "This method was mis-named, use SetMaxAng() instead!");
1094 SetMaxAng(x);
1095}
1096////////////////////////////////////////////////////////////////////////////////
1097/// Get maximum step angle.
1098/// WARNING -- this method / variable was mis-named.
1099
1101{
1102 Warning("GetMinAng", "This method was mis-named, use GetMaxAng() instead!");
1103 return GetMaxAng();
1104}
1105
1106////////////////////////////////////////////////////////////////////////////////
1107/// Set maximum step angle and rebuild tracks.
1108
1110{
1111 fH.fMaxAng = x;
1112 RebuildTracks();
1113}
1114
1115////////////////////////////////////////////////////////////////////////////////
1116/// Set maximum step-size and rebuild tracks.
1117
1119{
1120 fH.fMaxStep = x;
1121 RebuildTracks();
1122}
1123
1124////////////////////////////////////////////////////////////////////////////////
1125/// Set maximum error and rebuild tracks.
1126
1128{
1129 fH.fDelta = x;
1130 RebuildTracks();
1131}
1132
1133////////////////////////////////////////////////////////////////////////////////
1134/// Set daughter creation point fitting and rebuild tracks.
1135
1137{
1138 fFitDaughters = x;
1139 RebuildTracks();
1140}
1141
1142////////////////////////////////////////////////////////////////////////////////
1143/// Set track-reference fitting and rebuild tracks.
1144
1146{
1147 fFitReferences = x;
1148 RebuildTracks();
1149}
1150
1151////////////////////////////////////////////////////////////////////////////////
1152/// Set decay fitting and rebuild tracks.
1153
1155{
1156 fFitDecay = x;
1157 RebuildTracks();
1158}
1159
1160////////////////////////////////////////////////////////////////////////////////
1161/// Set line segment fitting and rebuild tracks.
1162
1164{
1166 RebuildTracks();
1167}
1168
1169////////////////////////////////////////////////////////////////////////////////
1170/// Set 2D-cluster fitting and rebuild tracks.
1171
1173{
1174 fFitCluster2Ds = x;
1175 RebuildTracks();
1176}
1177
1178////////////////////////////////////////////////////////////////////////////////
1179/// Set decay rendering and rebuild tracks.
1180
1182{
1183 fRnrDecay = rnr;
1184 RebuildTracks();
1185}
1186
1187////////////////////////////////////////////////////////////////////////////////
1188/// Set rendering of 2D-clusters and rebuild tracks.
1189
1191{
1192 fRnrCluster2Ds = rnr;
1193 RebuildTracks();
1194}
1195
1196////////////////////////////////////////////////////////////////////////////////
1197/// Set daughter rendering and rebuild tracks.
1198
1200{
1201 fRnrDaughters = rnr;
1202 RebuildTracks();
1203}
1204
1205////////////////////////////////////////////////////////////////////////////////
1206/// Set track-reference rendering and rebuild tracks.
1207
1209{
1210 fRnrReferences = rnr;
1211 RebuildTracks();
1212}
1213
1214////////////////////////////////////////////////////////////////////////////////
1215/// Set first-vertex rendering and rebuild tracks.
1216
1218{
1219 fRnrFV = x;
1220 RebuildTracks();
1221}
1222
1223////////////////////////////////////////////////////////////////////////////////
1224/// Set projection break-point mode and rebuild tracks.
1225
1227{
1229 RebuildTracks();
1230}
1231
1232////////////////////////////////////////////////////////////////////////////////
1233/// Set projection break-point rendering and rebuild tracks.
1234
1236{
1237 fRnrPTBMarkers = x;
1238 RebuildTracks();
1239}
1240
1241////////////////////////////////////////////////////////////////////////////////
1242/// Wrapper to step with method RungeKutta.
1243
1245 Double_t* vect, Double_t* vout)
1246{
1247 /// ******************************************************************
1248 /// * *
1249 /// * Runge-Kutta method for tracking a particle through a magnetic *
1250 /// * field. Uses Nystroem algorithm (See Handbook Nat. Bur. of *
1251 /// * Standards, procedure 25.5.20) *
1252 /// * *
1253 /// * Input parameters *
1254 /// * CHARGE Particle charge *
1255 /// * STEP Step size *
1256 /// * VECT Initial co-ords,direction cosines,momentum *
1257 /// * Output parameters *
1258 /// * VOUT Output co-ords,direction cosines,momentum *
1259 /// * User routine called *
1260 /// * CALL GUFLD(X,F) *
1261 /// * *
1262 /// * ==>Called by : <USER>, GUSWIM *
1263 /// * Authors R.Brun, M.Hansroul ********* *
1264 /// * V.Perevoztchikov (CUT STEP implementation) *
1265 /// * *
1266 /// * *
1267 /// ******************************************************************
1268
1269 Double_t h2, h4, f[4];
1270 Double_t /* xyzt[3], */ a, b, c, ph,ph2;
1271 Double_t secxs[4],secys[4],seczs[4],hxp[3];
1272 Double_t g1, g2, g3, g4, g5, g6, ang2, dxt, dyt, dzt;
1273 Double_t est, at, bt, ct, cba;
1274 Double_t f1, f2, f3, f4, rho, tet, hnorm, hp, rho1, sint, cost;
1275
1276 Double_t x;
1277 Double_t y;
1278 Double_t z;
1279
1280 Double_t xt;
1281 Double_t yt;
1282 Double_t zt;
1283
1284 // const Int_t maxit = 1992;
1285 const Int_t maxit = 500;
1286 const Int_t maxcut = 11;
1287
1288 const Double_t hmin = 1e-4; // !!! MT ADD, should be member
1289 const Double_t kdlt = 1e-3; // !!! MT CHANGE from 1e-4, should be member
1290 const Double_t kdlt32 = kdlt/32.;
1291 const Double_t kthird = 1./3.;
1292 const Double_t khalf = 0.5;
1293 const Double_t kec = 2.9979251e-3;
1294
1295 const Double_t kpisqua = 9.86960440109;
1296 const Int_t kix = 0;
1297 const Int_t kiy = 1;
1298 const Int_t kiz = 2;
1299 const Int_t kipx = 3;
1300 const Int_t kipy = 4;
1301 const Int_t kipz = 5;
1302
1303 // *.
1304 // *. ------------------------------------------------------------------
1305 // *.
1306 // * this constant is for units cm,gev/c and kgauss
1307 // *
1308 Int_t iter = 0;
1309 Int_t ncut = 0;
1310 for(Int_t j = 0; j < 7; j++)
1311 vout[j] = vect[j];
1312
1313 Double_t pinv = kec * fH.fCharge / vect[6];
1314 Double_t tl = 0.;
1315 Double_t h = step;
1316 Double_t rest;
1317
1318 do {
1319 rest = step - tl;
1320 if (TMath::Abs(h) > TMath::Abs(rest))
1321 h = rest;
1322
1323 f[0] = -fH.fB.fX;
1324 f[1] = -fH.fB.fY;
1325 f[2] = -fH.fB.fZ;
1326
1327 // * start of integration
1328 x = vout[0];
1329 y = vout[1];
1330 z = vout[2];
1331 a = vout[3];
1332 b = vout[4];
1333 c = vout[5];
1334
1335 h2 = khalf * h;
1336 h4 = khalf * h2;
1337 ph = pinv * h;
1338 ph2 = khalf * ph;
1339 secxs[0] = (b * f[2] - c * f[1]) * ph2;
1340 secys[0] = (c * f[0] - a * f[2]) * ph2;
1341 seczs[0] = (a * f[1] - b * f[0]) * ph2;
1342 ang2 = (secxs[0]*secxs[0] + secys[0]*secys[0] + seczs[0]*seczs[0]);
1343 if (ang2 > kpisqua) break;
1344
1345 dxt = h2 * a + h4 * secxs[0];
1346 dyt = h2 * b + h4 * secys[0];
1347 dzt = h2 * c + h4 * seczs[0];
1348 xt = x + dxt;
1349 yt = y + dyt;
1350 zt = z + dzt;
1351
1352 // * second intermediate point
1353 est = TMath::Abs(dxt) + TMath::Abs(dyt) + TMath::Abs(dzt);
1354 if (est > h) {
1355 if (ncut++ > maxcut) break;
1356 h *= khalf;
1357 continue;
1358 }
1359
1360 // xyzt[0] = xt;
1361 // xyzt[1] = yt;
1362 // xyzt[2] = zt;
1363
1364 fH.fB = fMagFieldObj->GetFieldD(xt, yt, zt);
1365 f[0] = -fH.fB.fX;
1366 f[1] = -fH.fB.fY;
1367 f[2] = -fH.fB.fZ;
1368
1369 at = a + secxs[0];
1370 bt = b + secys[0];
1371 ct = c + seczs[0];
1372
1373 secxs[1] = (bt * f[2] - ct * f[1]) * ph2;
1374 secys[1] = (ct * f[0] - at * f[2]) * ph2;
1375 seczs[1] = (at * f[1] - bt * f[0]) * ph2;
1376 at = a + secxs[1];
1377 bt = b + secys[1];
1378 ct = c + seczs[1];
1379 secxs[2] = (bt * f[2] - ct * f[1]) * ph2;
1380 secys[2] = (ct * f[0] - at * f[2]) * ph2;
1381 seczs[2] = (at * f[1] - bt * f[0]) * ph2;
1382 dxt = h * (a + secxs[2]);
1383 dyt = h * (b + secys[2]);
1384 dzt = h * (c + seczs[2]);
1385 xt = x + dxt;
1386 yt = y + dyt;
1387 zt = z + dzt;
1388 at = a + 2.*secxs[2];
1389 bt = b + 2.*secys[2];
1390 ct = c + 2.*seczs[2];
1391
1392 est = TMath::Abs(dxt)+TMath::Abs(dyt)+TMath::Abs(dzt);
1393 if (est > 2.*TMath::Abs(h)) {
1394 if (ncut++ > maxcut) break;
1395 h *= khalf;
1396 continue;
1397 }
1398
1399 // xyzt[0] = xt;
1400 // xyzt[1] = yt;
1401 // xyzt[2] = zt;
1402
1403 fH.fB = fMagFieldObj->GetFieldD(xt, yt, zt);
1404 f[0] = -fH.fB.fX;
1405 f[1] = -fH.fB.fY;
1406 f[2] = -fH.fB.fZ;
1407
1408 z = z + (c + (seczs[0] + seczs[1] + seczs[2]) * kthird) * h;
1409 y = y + (b + (secys[0] + secys[1] + secys[2]) * kthird) * h;
1410 x = x + (a + (secxs[0] + secxs[1] + secxs[2]) * kthird) * h;
1411
1412 secxs[3] = (bt*f[2] - ct*f[1])* ph2;
1413 secys[3] = (ct*f[0] - at*f[2])* ph2;
1414 seczs[3] = (at*f[1] - bt*f[0])* ph2;
1415 a = a+(secxs[0]+secxs[3]+2. * (secxs[1]+secxs[2])) * kthird;
1416 b = b+(secys[0]+secys[3]+2. * (secys[1]+secys[2])) * kthird;
1417 c = c+(seczs[0]+seczs[3]+2. * (seczs[1]+seczs[2])) * kthird;
1418
1419 est = TMath::Abs(secxs[0]+secxs[3] - (secxs[1]+secxs[2]))
1420 + TMath::Abs(secys[0]+secys[3] - (secys[1]+secys[2]))
1421 + TMath::Abs(seczs[0]+seczs[3] - (seczs[1]+seczs[2]));
1422
1423 if (est > kdlt && TMath::Abs(h) > hmin) {
1424 if (ncut++ > maxcut) break;
1425 h *= khalf;
1426 continue;
1427 }
1428
1429 ncut = 0;
1430 // * if too many iterations, go to helix
1431 if (iter++ > maxit) break;
1432
1433 tl += h;
1434 if (est < kdlt32)
1435 h *= 2.;
1436 cba = 1./ TMath::Sqrt(a*a + b*b + c*c);
1437 vout[0] = x;
1438 vout[1] = y;
1439 vout[2] = z;
1440 vout[3] = cba*a;
1441 vout[4] = cba*b;
1442 vout[5] = cba*c;
1443 rest = step - tl;
1444 if (step < 0.) rest = -rest;
1445 if (rest < 1.e-5*TMath::Abs(step))
1446 {
1447 Float_t dot = (vout[3]*vect[3] + vout[4]*vect[4] + vout[5]*vect[5]);
1448 fH.fPhi += TMath::ACos(dot);
1449 return;
1450 }
1451
1452 } while(1);
1453
1454 // angle too big, use helix
1455
1456 f1 = f[0];
1457 f2 = f[1];
1458 f3 = f[2];
1459 f4 = TMath::Sqrt(f1*f1+f2*f2+f3*f3);
1460 rho = -f4*pinv;
1461 tet = rho * step;
1462
1463 hnorm = 1./f4;
1464 f1 = f1*hnorm;
1465 f2 = f2*hnorm;
1466 f3 = f3*hnorm;
1467
1468 hxp[0] = f2*vect[kipz] - f3*vect[kipy];
1469 hxp[1] = f3*vect[kipx] - f1*vect[kipz];
1470 hxp[2] = f1*vect[kipy] - f2*vect[kipx];
1471
1472 hp = f1*vect[kipx] + f2*vect[kipy] + f3*vect[kipz];
1473
1474 rho1 = 1./rho;
1475 sint = TMath::Sin(tet);
1476 cost = 2.*TMath::Sin(khalf*tet)*TMath::Sin(khalf*tet);
1477
1478 g1 = sint*rho1;
1479 g2 = cost*rho1;
1480 g3 = (tet-sint) * hp*rho1;
1481 g4 = -cost;
1482 g5 = sint;
1483 g6 = cost * hp;
1484
1485 vout[kix] = vect[kix] + g1*vect[kipx] + g2*hxp[0] + g3*f1;
1486 vout[kiy] = vect[kiy] + g1*vect[kipy] + g2*hxp[1] + g3*f2;
1487 vout[kiz] = vect[kiz] + g1*vect[kipz] + g2*hxp[2] + g3*f3;
1488
1489 vout[kipx] = vect[kipx] + g4*vect[kipx] + g5*hxp[0] + g6*f1;
1490 vout[kipy] = vect[kipy] + g4*vect[kipy] + g5*hxp[1] + g6*f2;
1491 vout[kipz] = vect[kipz] + g4*vect[kipz] + g5*hxp[2] + g6*f3;
1492
1493 fH.fPhi += tet;
1494}
ROOT::R::TRInterface & r
Definition: Object.C:4
#define d(i)
Definition: RSha256.hxx:102
#define b(i)
Definition: RSha256.hxx:100
#define f(i)
Definition: RSha256.hxx:104
#define c(i)
Definition: RSha256.hxx:101
#define h(i)
Definition: RSha256.hxx:106
#define e(i)
Definition: RSha256.hxx:103
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
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
float Float_t
Definition: RtypesCore.h:57
const Bool_t kTRUE
Definition: RtypesCore.h:100
#define ClassImp(name)
Definition: Rtypes.h:364
@ kRed
Definition: Rtypes.h:66
@ kBlue
Definition: Rtypes.h:66
@ kYellow
Definition: Rtypes.h:66
Double_t Dot(const TGLVector3 &v1, const TGLVector3 &v2)
Definition: TGLUtil.h:317
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
Definition: TAttMarker.h:38
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition: TAttMarker.h:40
virtual void SetMarkerSize(Size_t msize=1)
Set the marker size.
Definition: TAttMarker.h:41
A list of TEveElements.
Definition: TEveElement.h:433
void StampObjProps()
Definition: TEveElement.h:399
virtual void CheckReferenceCount(const TEveException &eh="TEveElement::CheckReferenceCount ")
Check external references to this and eventually auto-destruct the render-element.
virtual void ElementChanged(Bool_t update_scenes=kTRUE, Bool_t redraw=kFALSE)
Call this after an element has been changed so that the state can be propagated around the framework.
Exception class thrown by TEve classes and macros.
Definition: TEveUtil.h:102
Implements constant magnetic field, given by a vector fB.
Implements constant magnetic filed that switches on given axial radius fR2 from vector fBIn to fBOut.
Abstract base-class for interfacing to magnetic field needed by the TEveTrackPropagator.
virtual Bool_t IsConst() const
virtual void PrintField(Double_t x, Double_t y, Double_t z) const
virtual TEveVectorD GetFieldD(Double_t x, Double_t y, Double_t z) const
virtual Double_t GetMaxFieldMagD() const
TEvePointSet is a render-element holding a collection of 3D points with optional per-point TRef and a...
Definition: TEvePointSet.h:36
Base-class for reference-counted objects with reverse references to TEveElement objects.
Definition: TEveUtil.h:187
RefMap_t::iterator RefMap_i
Definition: TEveUtil.h:190
RefMap_t fBackRefs
Definition: TEveUtil.h:192
Int_t fRefCount
Definition: TEveUtil.h:165
Holding structure for a number of track rendering parameters.
Double_t GetMaxAng() const
virtual void GoToBounds(TEveVectorD &p)
Propagate particle to bounds.
void SetFitReferences(Bool_t x)
Set track-reference fitting and rebuild tracks.
void SetRnrDecay(Bool_t x)
Set decay rendering and rebuild tracks.
void SetRnrDaughters(Bool_t x)
Set daughter rendering and rebuild tracks.
void SetFitLineSegments(Bool_t x)
Set line segment fitting and rebuild tracks.
void ResetTrack()
Reset cache holding particle trajectory.
void SetMagFieldObj(TEveMagField *field, Bool_t own_field=kTRUE)
Set constant magnetic field and rebuild tracks.
Bool_t HelixIntersectPlane(const TEveVectorD &p, const TEveVectorD &point, const TEveVectorD &normal, TEveVectorD &itsect)
Intersect helix with a plane.
void StepRungeKutta(Double_t step, Double_t *vect, Double_t *vout)
Wrapper to step with method RungeKutta.
void SetDelta(Double_t x)
Set maximum error and rebuild tracks.
Bool_t IntersectPlane(const TEveVectorD &p, const TEveVectorD &point, const TEveVectorD &normal, TEveVectorD &itsect)
Find intersection of currently propagated track with a plane.
void DistributeOffset(const TEveVectorD &off, Int_t first_point, Int_t np, TEveVectorD &p)
Distribute offset between first and last point index and rotate momentum.
static Bool_t IsOutsideBounds(const TEveVectorD &point, Double_t maxRsqr, Double_t maxZ)
void ClosestPointFromVertexToLineSegment(const TEveVectorD &v, const TEveVectorD &s, const TEveVectorD &r, Double_t rMagInv, TEveVectorD &c)
Get closest point from given vertex v to line segment defined with s and r.
void Step(const TEveVector4D &v, const TEveVectorD &p, TEveVector4D &vOut, TEveVectorD &pOut)
Wrapper to step helix.
void SetMaxR(Double_t x)
Set maximum radius and rebuild tracks.
Double_t GetTrackLength(Int_t start_point=0, Int_t end_point=-1) const
Calculate track length from start_point to end_point.
void SetFitDaughters(Bool_t x)
Set daughter creation point fitting and rebuild tracks.
TEveTrackPropagator(const TEveTrackPropagator &)
virtual Bool_t GoToVertex(TEveVectorD &v, TEveVectorD &p)
Propagate particle with momentum p to vertex v.
static Double_t fgEditorMaxZ
Bool_t PointOverVertex(const TEveVector4D &v0, const TEveVector4D &v, Double_t *p=0)
void LineToBounds(TEveVectorD &p)
Propagate neutral particle with momentum p to bounds.
Int_t GetCurrentPoint() const
Get index of current point on track.
virtual ~TEveTrackPropagator()
Destructor.
static Double_t fgDefMagField
void Update(const TEveVector4D &v, const TEveVectorD &p, Bool_t full_update=kFALSE, Bool_t enforce_max_step=kFALSE)
Update helix / B-field projection state.
void FillPointSet(TEvePointSet *ps) const
Reset ps and populate it with points in propagation cache.
virtual void OnZeroRefCount()
Virtual from TEveRefBackPtr - track reference count has reached zero.
void SetRnrFV(Bool_t x)
Set first-vertex rendering and rebuild tracks.
void SetProjTrackBreaking(UChar_t x)
Set projection break-point mode and rebuild tracks.
Bool_t ClosestPointBetweenLines(const TEveVectorD &, const TEveVectorD &, const TEveVectorD &, const TEveVectorD &, TEveVectorD &out)
Get closest point on line defined with vector p0 and u.
void SetMinAng(Double_t x)
Set maximum step angle and rebuild tracks.
void LoopToBounds(TEveVectorD &p)
Propagate charged particle with momentum p to bounds.
static TEveTrackPropagator fgDefault
Bool_t LineIntersectPlane(const TEveVectorD &p, const TEveVectorD &point, const TEveVectorD &normal, TEveVectorD &itsect)
Intersect line with a plane.
TEveMagField * fMagFieldObj
void SetRnrCluster2Ds(Bool_t x)
Set rendering of 2D-clusters and rebuild tracks.
std::vector< TEveVector4D > fPoints
virtual Bool_t GoToLineSegment(const TEveVectorD &s, const TEveVectorD &r, TEveVectorD &p)
Propagate particle with momentum p to line with start point s and vector r to the second point.
void SetRnrReferences(Bool_t x)
Set track-reference rendering and rebuild tracks.
static Double_t fgEditorMaxR
void SetMaxAng(Double_t x)
Set maximum step angle and rebuild tracks.
Bool_t LineToVertex(TEveVectorD &v)
Propagate neutral particle to vertex v.
void InitTrack(const TEveVectorD &v, Int_t charge)
Initialize internal data-members for given particle parameters.
virtual void ElementChanged(Bool_t update_scenes=kTRUE, Bool_t redraw=kFALSE)
Element-change notification.
static const Double_t fgkB2C
void SetMaxStep(Double_t x)
Set maximum step-size and rebuild tracks.
void SetFitCluster2Ds(Bool_t x)
Set 2D-cluster fitting and rebuild tracks.
std::vector< TEveVector4D > fLastPoints
Bool_t LoopToLineSegment(const TEveVectorD &s, const TEveVectorD &r, TEveVectorD &p)
Propagate charged particle with momentum p to line segment with point s and vector r to the second po...
Double_t GetMinAng() const
Get maximum step angle.
void PrintMagField(Double_t x, Double_t y, Double_t z) const
void RebuildTracks()
Rebuild all tracks using this render-style.
void SetRnrPTBMarkers(Bool_t x)
Set projection break-point rendering and rebuild tracks.
virtual void CheckReferenceCount(const TEveException &eh="TEveElement::CheckReferenceCount ")
Check reference count - virtual from TEveElement.
Bool_t LoopToVertex(TEveVectorD &v, TEveVectorD &p)
Propagate charged particle with momentum p to vertex v.
void SetMaxZ(Double_t x)
Set maximum z and rebuild tracks.
void SetFitDecay(Bool_t x)
Set decay fitting and rebuild tracks.
void SetMaxOrbs(Double_t x)
Set maximum number of orbits and rebuild tracks.
void SetMagField(Double_t bX, Double_t bY, Double_t bZ)
Set constant magnetic field and rebuild tracks.
Visual representation of a track.
Definition: TEveTrack.h:33
virtual void MakeTrack(Bool_t recurse=kTRUE)
Calculate track representation based on track data and current settings of the propagator.
Definition: TEveTrack.cxx:336
TEveTrans is a 4x4 transformation matrix for homogeneous coordinates stored internally in a column-ma...
Definition: TEveTrans.h:27
Minimal, templated four-vector.
Definition: TEveVector.h:243
TT Perp2() const
Definition: TEveVector.h:101
TT Normalize(TT length=1)
Normalize the vector to length if current length is non-zero.
Definition: TEveVector.cxx:56
TT Dot(const TEveVectorT &a) const
Definition: TEveVector.h:168
TT Mag() const
Definition: TEveVector.h:99
TT Mag2() const
Definition: TEveVector.h:98
TT R() const
Definition: TEveVector.h:103
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:879
T Mag(const SVector< T, D > &rhs)
Vector magnitude (Euclidian norm) Compute : .
Definition: Functions.h:253
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
TF1 * f1
Definition: legend1.C:11
Double_t Sqrt(Double_t x)
static constexpr double nm
static constexpr double s
static constexpr double ps
TMath.
Definition: TMathBase.h:35
Double_t ACos(Double_t)
Definition: TMath.h:669
Short_t Range(Short_t lb, Short_t ub, Short_t x)
Definition: TMathBase.h:240
constexpr Double_t DegToRad()
Conversion from degree to radian:
Definition: TMath.h:81
Double_t Sqrt(Double_t x)
Definition: TMath.h:691
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:176
Double_t Cos(Double_t)
Definition: TMath.h:643
Double_t Sin(Double_t)
Definition: TMath.h:639
T * Cross(const T v1[3], const T v2[3], T out[3])
Calculate the Cross Product of two vectors: out = [v1 x v2].
Definition: TMath.h:1175
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
constexpr Double_t TwoPi()
Definition: TMath.h:44
void UpdateRK(const TEveVectorD &p, const TEveVectorD &b)
Update helix for stepper RungeKutta.
void UpdateHelix(const TEveVectorD &p, const TEveVectorD &b, Bool_t full_update, Bool_t enforce_max_step)
Update helix parameters.
void Step(const TEveVector4D &v, const TEveVectorD &p, TEveVector4D &vOut, TEveVectorD &pOut)
Step helix for given momentum p from vertex v.
void UpdateCommon(const TEveVectorD &p, const TEveVectorD &b)
Common update code for helix and RK propagation.
auto * tt
Definition: textangle.C:16
auto * a
Definition: textangle.C:12
static uint64_t sum(uint64_t i)
Definition: Factory.cxx:2345