Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
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
39
40/** \class TEveMagFieldConst
41\ingroup TEve
42Implements constant magnetic field, given by a vector fB.
43
44NOTE: Magnetic field direction convention is inverted.
45*/
46
47
48/** \class TEveMagFieldDuo
49\ingroup TEve
50Implements constant magnetic filed that switches on given axial radius fR2
51from vector fBIn to fBOut.
52
53NOTE: Magnetic field direction convention is inverted.
54*/
55
56
57namespace
58{
59 //const Double_t kBMin = 1e-6;
60 const Double_t kPtMinSqr = 1e-20;
61 const Double_t kAMin = 1e-10;
62 const Double_t kStepEps = 1e-3;
63}
64
65////////////////////////////////////////////////////////////////////////////////
66/// Default constructor.
67
69 fCharge(0),
70 fMaxAng(45), fMaxStep(20.f), fDelta(0.1),
71 fPhi(0), fValid(kFALSE),
72 fLam(-1), fR(-1), fPhiStep(-1), fSin(-1), fCos(-1),
73 fRKStep(20.0),
74 fPtMag(-1), fPlMag(-1), fLStep(-1)
75{
76}
77
78////////////////////////////////////////////////////////////////////////////////
79/// Common update code for helix and RK propagation.
80
82{
83 fB = b;
84
85 // base vectors
86 fE1 = b;
87 fE1.Normalize();
88 fPlMag = p.Dot(fE1);
89 fPl = fE1*fPlMag;
90
91 fPt = p - fPl;
92 fPtMag = fPt.Mag();
93 fE2 = fPt;
94 fE2.Normalize();
95}
96
97////////////////////////////////////////////////////////////////////////////////
98/// Update helix parameters.
99
102{
103 UpdateCommon(p, b);
104
105 // helix parameters
106 TMath::Cross(fE1.Arr(), fE2.Arr(), fE3.Arr());
107 if (fCharge < 0) fE3.NegateXYZ();
108
109 if (full_update)
110 {
111 using namespace TMath;
112 Double_t a = fgkB2C * b.Mag() * Abs(fCharge);
113 if (a > kAMin && fPtMag*fPtMag > kPtMinSqr)
114 {
115 fValid = kTRUE;
116
117 fR = Abs(fPtMag / a);
118 fLam = fPlMag / fPtMag;
119
120 // get phi step, compare fMaxAng with fDelta
121 fPhiStep = fMaxAng * DegToRad();
122 if (fR > fDelta)
123 {
124 Double_t ang = 2.0 * ACos(1.0f - fDelta/fR);
125 if (ang < fPhiStep)
126 fPhiStep = ang;
127 }
128
129 // check max step size
130 Double_t curr_step = fR * fPhiStep * Sqrt(1.0f + fLam*fLam);
131 if (curr_step > fMaxStep || enforce_max_step)
132 fPhiStep *= fMaxStep / curr_step;
133
134 fLStep = fR * fPhiStep * fLam;
135 fSin = Sin(fPhiStep);
136 fCos = Cos(fPhiStep);
137 }
138 else
139 {
140 fValid = kFALSE;
141 }
142 }
143}
144
145////////////////////////////////////////////////////////////////////////////////
146/// Update helix for stepper RungeKutta.
147
149{
150 UpdateCommon(p, b);
151
152 if (fCharge)
153 {
154 fValid = kTRUE;
155
156 // cached values for propagator
157 fB = b;
158 fPlMag = p.Dot(fB);
159 }
160 else
161 {
162 fValid = kFALSE;
163 }
164}
165
166////////////////////////////////////////////////////////////////////////////////
167/// Step helix for given momentum p from vertex v.
168
171{
172 vOut = v;
173
174 if (fValid)
175 {
176 TEveVectorD d = fE2*(fR*fSin) + fE3*(fR*(1-fCos)) + fE1*fLStep;
177 vOut += d;
178 vOut.fT += TMath::Abs(fLStep);
179
180 pOut = fPl + fE2*(fPtMag*fCos) + fE3*(fPtMag*fSin);
181
182 fPhi += fPhiStep;
183 }
184 else
185 {
186 // case: pT < kPtMinSqr or B < kBMin
187 // might happen if field directon changes pT ~ 0 or B becomes zero
188 vOut += p * (fMaxStep / p.Mag());
189 vOut.fT += fMaxStep;
190 pOut = p;
191 }
192}
193
194/** \class TEveTrackPropagator
195\ingroup TEve
196Holding structure for a number of track rendering parameters.
197Calculates path taking into account the parameters.
198
199NOTE: Magnetic field direction convention is inverted.
200
201This is decoupled from TEveTrack/TEveTrackList to allow sharing of the
202Propagator among several instances. Back references are kept so the tracks
203can be recreated when the parameters change.
204
205TEveTrackList has Get/Set methods for RnrStlye. TEveTrackEditor and
206TEveTrackListEditor provide editor access.
207
208Enum EProjTrackBreaking_e and member fProjTrackBreaking specify whether 2D
209projected tracks get broken into several segments when the projected space
210consists of separate domains (like Rho-Z). The track-breaking is enabled by
211default.
212*/
213
214
216const Double_t TEveTrackPropagator::fgkB2C = 0.299792458e-2;
218
221
222////////////////////////////////////////////////////////////////////////////////
223/// Default constructor.
224
267
268////////////////////////////////////////////////////////////////////////////////
269/// Destructor.
270
272{
273 if (fOwnMagFiledObj)
274 {
275 delete fMagFieldObj;
276 }
277}
278
279////////////////////////////////////////////////////////////////////////////////
280/// Virtual from TEveRefBackPtr - track reference count has reached zero.
281
283{
284 CheckReferenceCount("TEveTrackPropagator::OnZeroRefCount ");
285}
286
287////////////////////////////////////////////////////////////////////////////////
288/// Check reference count - virtual from TEveElement.
289/// Must also take into account references from TEveRefBackPtr.
290
298
299////////////////////////////////////////////////////////////////////////////////
300/// Element-change notification.
301/// Stamp all tracks as requiring display-list regeneration.
302/// Virtual from TEveElement.
303
305{
307 RefMap_i i = fBackRefs.begin();
308 while (i != fBackRefs.end())
309 {
310 track = dynamic_cast<TEveTrack*>(i->first);
311 track->StampObjProps();
312 ++i;
313 }
315}
316
317////////////////////////////////////////////////////////////////////////////////
318/// Initialize internal data-members for given particle parameters.
319
321{
322 fV = v;
323 fPoints.push_back(fV);
324
325 // init helix
326 fH.fPhi = 0;
327 fH.fCharge = charge;
328}
329
330////////////////////////////////////////////////////////////////////////////////
331/// TEveVectorF wrapper.
332
338
339////////////////////////////////////////////////////////////////////////////////
340/// Reset cache holding particle trajectory.
341
343{
344 fLastPoints.clear();
345 fPoints.swap(fLastPoints);
346
347 // reset helix
348 fH.fPhi = 0;
349}
350
351////////////////////////////////////////////////////////////////////////////////
352/// Get index of current point on track.
353
355{
356 return fPoints.size() - 1;
357}
358
359////////////////////////////////////////////////////////////////////////////////
360/// Calculate track length from start_point to end_point.
361/// If end_point is less than 0, distance to the end is returned.
362
364{
365 if (end_point < 0) end_point = fPoints.size() - 1;
366
367 Double_t sum = 0;
368 for (Int_t i = start_point; i < end_point; ++i)
369 {
370 sum += (fPoints[i+1] - fPoints[i]).Mag();
371 }
372 return sum;
373}
374
375////////////////////////////////////////////////////////////////////////////////
376/// Propagate particle with momentum p to vertex v.
377
379{
380 Update(fV, p, kTRUE);
381
382 if ((v-fV).Mag() < kStepEps)
383 {
384 fPoints.push_back(v);
385 return kTRUE;
386 }
387
388 return fH.fValid ? LoopToVertex(v, p) : LineToVertex(v);
389}
390
391////////////////////////////////////////////////////////////////////////////////
392/// Propagate particle with momentum p to line with start point s and vector r
393/// to the second point.
394
396{
397 Update(fV, p, kTRUE);
398
399 if (!fH.fValid)
400 {
404 return kTRUE;
405 }
406 else
407 {
408 return LoopToLineSegment(s, r, p);
409 }
410}
411
412////////////////////////////////////////////////////////////////////////////////
413/// TEveVectorF wrapper.
414
422
423////////////////////////////////////////////////////////////////////////////////
424/// TEveVectorF wrapper.
425
433
434////////////////////////////////////////////////////////////////////////////////
435/// Propagate particle to bounds.
436/// Return TRUE if hit bounds.
437
444
445////////////////////////////////////////////////////////////////////////////////
446/// TEveVectorF wrapper.
447
454
455////////////////////////////////////////////////////////////////////////////////
456/// Update helix / B-field projection state.
457
460{
461 if (fStepper == kHelix)
462 {
464 }
465 else
466 {
468
469 if (full_update)
470 {
471 using namespace TMath;
472
474 if (a > kAMin)
475 {
476 fH.fR = p.Mag() / a;
477
478 // get phi step, compare fDelta with MaxAng
479 fH.fPhiStep = fH.fMaxAng * DegToRad();
480 if (fH.fR > fH.fDelta )
481 {
482 Double_t ang = 2.0 * ACos(1.0f - fH.fDelta/fH.fR);
483 if (ang < fH.fPhiStep)
484 fH.fPhiStep = ang;
485 }
486
487 // check against maximum step-size
488 fH.fRKStep = fH.fR * fH.fPhiStep * Sqrt(1 + fH.fLam*fH.fLam);
490 {
493 }
494 }
495 else
496 {
498 }
499 }
500 }
501}
502
503////////////////////////////////////////////////////////////////////////////////
504/// Wrapper to step helix.
505
507{
508 if (fStepper == kHelix)
509 {
510 fH.Step(v, p, vOut, pOut);
511 }
512 else
513 {
514 Double_t vecRKIn[7];
515 vecRKIn[0] = v.fX;
516 vecRKIn[1] = v.fY;
517 vecRKIn[2] = v.fZ;
518 Double_t pm = p.Mag();
519 Double_t nm = 1.0 / pm;
520 vecRKIn[3] = p.fX*nm;
521 vecRKIn[4] = p.fY*nm;
522 vecRKIn[5] = p.fZ*nm;
523 vecRKIn[6] = p.Mag();
524
527
528 vOut.fX = vecRKOut[0];
529 vOut.fY = vecRKOut[1];
530 vOut.fZ = vecRKOut[2];
531 vOut.fT = v.fT + fH.fRKStep;
532 pm = vecRKOut[6];
533 pOut.fX = vecRKOut[3]*pm;
534 pOut.fY = vecRKOut[4]*pm;
535 pOut.fZ = vecRKOut[5]*pm;
536 }
537}
538
539////////////////////////////////////////////////////////////////////////////////
540/// Propagate charged particle with momentum p to bounds.
541/// It is expected that Update() with full-update was called before.
542
544{
545 const Double_t maxRsq = fMaxR*fMaxR;
546
550
551 Int_t np = fPoints.size();
553
554 while (fH.fPhi < maxPhi && np<fNMax)
555 {
556 Step(currV, p, forwV, forwP);
557
558 // cross R
559 if (forwV.Perp2() > maxRsq)
560 {
561 Float_t t = (fMaxR - currV.R()) / (forwV.R() - currV.R());
562 if (t < 0 || t > 1)
563 {
564 Warning("HelixToBounds", "In MaxR crossing expected t>=0 && t<=1: t=%f, r1=%f, r2=%f, MaxR=%f.",
565 t, currV.R(), forwV.R(), fMaxR);
566 return;
567 }
569 d -= currV;
570 d *= t;
571 d += currV;
572 fPoints.push_back(d);
573 return;
574 }
575
576 // cross Z
577 else if (TMath::Abs(forwV.fZ) > fMaxZ)
578 {
579 Double_t t = (fMaxZ - TMath::Abs(currV.fZ)) / TMath::Abs((forwV.fZ - currV.fZ));
580 if (t < 0 || t > 1)
581 {
582 Warning("HelixToBounds", "In MaxZ crossing expected t>=0 && t<=1: t=%f, z1=%f, z2=%f, MaxZ=%f.",
583 t, currV.fZ, forwV.fZ, fMaxZ);
584 return;
585 }
587 d *= t;
588 d += currV;
589 fPoints.push_back(d);
590 return;
591 }
592
593 currV = forwV;
594 p = forwP;
595 Update(currV, p);
596
597 fPoints.push_back(currV);
598 ++np;
599 }
600}
601
602////////////////////////////////////////////////////////////////////////////////
603/// Propagate charged particle with momentum p to vertex v.
604/// It is expected that Update() with full-update was called before.
605
607{
608 const Double_t maxRsq = fMaxR * fMaxR;
609
613
614 Int_t first_point = fPoints.size();
616
618
619 do
620 {
621 Step(currV, p, forwV, forwP);
623
625 {
626 break;
627 }
628
630 {
631 fV = currV;
632 return kFALSE;
633 }
634
635 fPoints.push_back(forwV);
636 currV = forwV;
637 p = forwP;
638 prod0 = prod1;
639 ++np;
640 } while (np < fNMax);
641
642 // make the remaining fractional step
643 if (np > first_point)
644 {
645 if ((v - currV).Mag() > kStepEps)
646 {
648 if (step_frac > 0)
649 {
650 // Step for fraction of previous step size.
651 // We pass 'enforce_max_step' flag to Update().
653 fH.fMaxStep = step_frac * (forwV - currV).Mag();
655 Step(currV, p, forwV, forwP);
656 p = forwP;
657 currV = forwV;
658 fPoints.push_back(currV);
659 ++np;
661 }
662
663 // Distribute offset to desired crossing point over all segment.
664
665 TEveVectorD off(v - currV);
666 off *= 1.0f / currV.fT;
668 fV = v;
669 return kTRUE;
670 }
671 }
672
673 fPoints.push_back(v);
674 fV = v;
675 return kTRUE;
676}
677
678////////////////////////////////////////////////////////////////////////////////
679/// Propagate charged particle with momentum p to line segment with point s and
680/// vector r to the second point. It is expected that Update() with full-update
681/// was called before. Returns kFALSE if hits bounds.
682
684{
685 const Double_t maxRsq = fMaxR * fMaxR;
686 const Double_t rMagInv = 1./r.Mag();
687
691
692 Int_t first_point = fPoints.size();
694
697 do
698 {
699 Step(currV, p, forwV, forwP);
701
703
704 // check forwV is over segment with orthogonal component of
705 // momentum to vector r
706 TEveVectorD b = r; b.Normalize();
707 Double_t x = forwP.Dot(b);
709 if (pTPM.Dot(forwC - forwV) < 0)
710 {
711 break;
712 }
713
715 {
716 fV = currV;
717 return kFALSE;
718 }
719
720 fPoints.push_back(forwV);
721 currV = forwV;
722 p = forwP;
723 currC = forwC;
724 ++np;
725 } while (np < fNMax);
726
727 // Get closest point on segment relative to line with forw and currV points.
730
731 // make the remaining fractional step
732 if (np > first_point)
733 {
734 if ((v - currV).Mag() > kStepEps)
735 {
737 TEveVector delta = v - currV;
738 Double_t step_frac = last_step.Dot(delta) / last_step.Mag2();
739 if (step_frac > 0)
740 {
741 // Step for fraction of previous step size.
742 // We pass 'enforce_max_step' flag to Update().
744 fH.fMaxStep = step_frac * (forwV - currV).Mag();
746 Step(currV, p, forwV, forwP);
747 p = forwP;
748 currV = forwV;
749 fPoints.push_back(currV);
750 ++np;
752 }
753
754 // Distribute offset to desired crossing point over all segment.
755
756 TEveVectorD off(v - currV);
757 off *= 1.0f / currV.fT;
759 fV = v;
760 return kTRUE;
761 }
762 }
763
764 fPoints.push_back(v);
765 fV = v;
766 return kTRUE;
767}
768
769////////////////////////////////////////////////////////////////////////////////
770/// Distribute offset between first and last point index and rotate
771/// momentum.
772
774{
775 // Calculate the required momentum rotation.
776 // lpd - last-points-delta
778 lpd0 -= fPoints[np-2];
779 lpd0.Normalize();
780
781 for (Int_t i = first_point; i < np; ++i)
782 {
783 fPoints[i] += off * fPoints[i].fT;
784 }
785
787 lpd1 -= fPoints[np-2];
788 lpd1.Normalize();
789
791 tt.SetupFromToVec(lpd0, lpd1);
792
793 // TEveVectorD pb4(p);
794 // printf("Rotating momentum: p0 = "); p.Dump();
795 tt.RotateIP(p);
796 // printf(" p1 = "); p.Dump();
797 // printf(" n1=%f, n2=%f, dp = %f deg\n", pb4.Mag(), p.Mag(),
798 // TMath::RadToDeg()*TMath::ACos(p.Dot(pb4)/(pb4.Mag()*p.Mag())));
799}
800
801////////////////////////////////////////////////////////////////////////////////
802/// Propagate neutral particle to vertex v.
803
805{
807
808 currV.fX = v.fX;
809 currV.fY = v.fY;
810 currV.fZ = v.fZ;
811 fPoints.push_back(currV);
812
813 fV = v;
814 return kTRUE;
815}
816
817////////////////////////////////////////////////////////////////////////////////
818/// Propagate neutral particle with momentum p to bounds.
819
821{
822 Double_t tZ = 0, tR = 0, tB = 0;
823
824 // time where particle intersect +/- fMaxZ
825 if (p.fZ > 0)
826 tZ = (fMaxZ - fV.fZ) / p.fZ;
827 else if (p.fZ < 0)
828 tZ = - (fMaxZ + fV.fZ) / p.fZ;
829 else
830 tZ = 1e99;
831
832 // time where particle intersects cylinder
833 Double_t a = p.fX*p.fX + p.fY*p.fY;
834 Double_t b = 2.0 * (fV.fX*p.fX + fV.fY*p.fY);
836 Double_t d = b*b - 4.0*a*c;
837 if (d >= 0) {
839 tR = (-b - sqrtD) / (2.0 * a);
840 if (tR < 0) {
841 tR = (-b + sqrtD) / (2.0 * a);
842 }
843 tB = tR < tZ ? tR : tZ; // compare the two times
844 } else {
845 tB = tZ;
846 }
847 TEveVectorD nv(fV.fX + p.fX*tB, fV.fY + p.fY*tB, fV.fZ + p.fZ*tB);
849}
850
851////////////////////////////////////////////////////////////////////////////////
852/// Intersect helix with a plane. Current position and argument p define
853/// the helix.
854
856 const TEveVectorD& point,
857 const TEveVectorD& normal,
859{
860 TEveVectorD pos(fV);
862 if (fMagFieldObj->IsConst())
864
866 TEveVectorD delta = pos - point;
867 Double_t d = delta.Dot(n);
868 if (d > 0) {
869 n.NegateXYZ(); // Turn normal around so that we approach from negative side of the plane
870 d = -d;
871 }
872
875 TEveVector4D pos4(pos);
876 while (kTRUE)
877 {
878 Update(pos4, mom);
879 Step(pos4, mom, forwV , forwP);
880 Double_t new_d = (forwV - point).Dot(n);
881 if (new_d < d)
882 {
883 // We are going further away ... fail intersect.
884 Warning("HelixIntersectPlane", "going away from the plane.");
885 return kFALSE;
886 }
887 if (new_d > 0)
888 {
889 delta = forwV - pos4;
890 itsect = pos4 + delta * ((point - pos4).Dot(n) / delta.Dot(n));
891
892 return kTRUE;
893 }
894 pos4 = forwV;
895 mom = forwP;
896 }
897}
898
899////////////////////////////////////////////////////////////////////////////////
900/// Intersect line with a plane. Current position and argument p define
901/// the line.
902
904 const TEveVectorD& point,
905 const TEveVectorD& normal,
907{
908 TEveVectorD pos(fV.fX, fV.fY, fV.fZ);
909 TEveVectorD delta = point - pos;
910
911 Double_t pn = p.Dot(normal);
912 if (pn == 0)
913 {
914 return kFALSE;
915 }
916 Double_t t = delta.Dot(normal) / pn;
917 if (t < 0) {
918 return kFALSE;
919 } else {
920 itsect = pos + p*t;
921 return kTRUE;
922 }
923}
924
925////////////////////////////////////////////////////////////////////////////////
926/// Find intersection of currently propagated track with a plane.
927/// Current track position is used as starting point.
928///
929/// Args:
930/// - p - track momentum to use for extrapolation
931/// - point - a point on a plane
932/// - normal - normal of the plane
933/// - itsect - output, point of intersection
934/// Returns:
935/// - kFALSE if intersection can not be found, kTRUE otherwise.
936
938 const TEveVectorD& point,
939 const TEveVectorD& normal,
941{
942 if (fH.fCharge && fMagFieldObj)
943 return HelixIntersectPlane(p, point, normal, itsect);
944 else
945 return LineIntersectPlane(p, point, normal, itsect);
946}
947
948////////////////////////////////////////////////////////////////////////////////
949/// Get closest point from given vertex v to line segment defined with s and r.
950/// Argument rMagInv is cached. rMagInv= 1./rMag()
951
953 const TEveVectorD& s,
954 const TEveVectorD& r,
956 TEveVectorD& c)
957{
958 TEveVectorD dir = v - s;
959 TEveVectorD b1 = r * rMagInv;
960
961 // parallel distance
962 Double_t dot = dir.Dot(b1);
963 TEveVectorD dirI = dot * b1;
964
965 Double_t facX = dot * rMagInv;
966
967 if (facX <= 0)
968 c = s;
969 else if (facX >= 1)
970 c = s + r;
971 else
972 c = s + dirI;
973}
974
975////////////////////////////////////////////////////////////////////////////////
976/// Get closest point on line defined with vector p0 and u.
977/// Return false if the point is forced on the line segment.
978
980 const TEveVectorD& u,
981 const TEveVectorD& q0,
982 const TEveVectorD& v,
983 TEveVectorD& out)
984{
985 TEveVectorD w0 = p0 -q0;
986 Double_t a = u.Mag2();
987 Double_t b = u.Dot(v);
988 Double_t c = v.Mag2();
989 Double_t d = u.Dot(w0);
990 Double_t e = v.Dot(w0);
991
992 Double_t x = (b*e - c*d)/(a*c -b*b);
994 out = p0 + TMath::Range(0., 1., x) * u;
995 return force;
996}
997
998////////////////////////////////////////////////////////////////////////////////
999/// Reset ps and populate it with points in propagation cache.
1000
1002{
1003 Int_t size = TMath::Min(fNMax, (Int_t)fPoints.size());
1004 ps->Reset(size);
1005 for (Int_t i = 0; i < size; ++i)
1006 {
1007 const TEveVector4D& v = fPoints[i];
1008 ps->SetNextPoint(v.fX, v.fY, v.fZ);
1009 }
1010}
1011
1012////////////////////////////////////////////////////////////////////////////////
1013/// Rebuild all tracks using this render-style.
1014
1016{
1018 RefMap_i i = fBackRefs.begin();
1019 while (i != fBackRefs.end())
1020 {
1021 track = dynamic_cast<TEveTrack*>(i->first);
1022 track->MakeTrack();
1023 track->StampObjProps();
1024 ++i;
1025 }
1026}
1027
1028////////////////////////////////////////////////////////////////////////////////
1029/// Set constant magnetic field and rebuild tracks.
1030
1035
1036////////////////////////////////////////////////////////////////////////////////
1037/// Set constant magnetic field and rebuild tracks.
1038
1048
1049////////////////////////////////////////////////////////////////////////////////
1050
1055
1056////////////////////////////////////////////////////////////////////////////////
1057/// Set maximum radius and rebuild tracks.
1058
1064
1065////////////////////////////////////////////////////////////////////////////////
1066/// Set maximum z and rebuild tracks.
1067
1073
1074////////////////////////////////////////////////////////////////////////////////
1075/// Set maximum number of orbits and rebuild tracks.
1076
1082
1083////////////////////////////////////////////////////////////////////////////////
1084/// Set maximum step angle and rebuild tracks.
1085/// WARNING -- this method / variable was mis-named.
1086
1088{
1089 Warning("SetMinAng", "This method was mis-named, use SetMaxAng() instead!");
1090 SetMaxAng(x);
1091}
1092////////////////////////////////////////////////////////////////////////////////
1093/// Get maximum step angle.
1094/// WARNING -- this method / variable was mis-named.
1095
1097{
1098 Warning("GetMinAng", "This method was mis-named, use GetMaxAng() instead!");
1099 return GetMaxAng();
1100}
1101
1102////////////////////////////////////////////////////////////////////////////////
1103/// Set maximum step angle and rebuild tracks.
1104
1110
1111////////////////////////////////////////////////////////////////////////////////
1112/// Set maximum step-size and rebuild tracks.
1113
1119
1120////////////////////////////////////////////////////////////////////////////////
1121/// Set maximum error and rebuild tracks.
1122
1128
1129////////////////////////////////////////////////////////////////////////////////
1130/// Set daughter creation point fitting and rebuild tracks.
1131
1137
1138////////////////////////////////////////////////////////////////////////////////
1139/// Set track-reference fitting and rebuild tracks.
1140
1146
1147////////////////////////////////////////////////////////////////////////////////
1148/// Set decay fitting and rebuild tracks.
1149
1155
1156////////////////////////////////////////////////////////////////////////////////
1157/// Set line segment fitting and rebuild tracks.
1158
1164
1165////////////////////////////////////////////////////////////////////////////////
1166/// Set 2D-cluster fitting and rebuild tracks.
1167
1173
1174////////////////////////////////////////////////////////////////////////////////
1175/// Set decay rendering and rebuild tracks.
1176
1182
1183////////////////////////////////////////////////////////////////////////////////
1184/// Set rendering of 2D-clusters and rebuild tracks.
1185
1191
1192////////////////////////////////////////////////////////////////////////////////
1193/// Set daughter rendering and rebuild tracks.
1194
1200
1201////////////////////////////////////////////////////////////////////////////////
1202/// Set track-reference rendering and rebuild tracks.
1203
1209
1210////////////////////////////////////////////////////////////////////////////////
1211/// Set first-vertex rendering and rebuild tracks.
1212
1218
1219////////////////////////////////////////////////////////////////////////////////
1220/// Set projection break-point mode and rebuild tracks.
1221
1227
1228////////////////////////////////////////////////////////////////////////////////
1229/// Set projection break-point rendering and rebuild tracks.
1230
1236
1237////////////////////////////////////////////////////////////////////////////////
1238/// Wrapper to step with method RungeKutta.
1239
1242{
1243 /// ******************************************************************
1244 /// * *
1245 /// * Runge-Kutta method for tracking a particle through a magnetic *
1246 /// * field. Uses Nystroem algorithm (See Handbook Nat. Bur. of *
1247 /// * Standards, procedure 25.5.20) *
1248 /// * *
1249 /// * Input parameters *
1250 /// * CHARGE Particle charge *
1251 /// * STEP Step size *
1252 /// * VECT Initial co-ords,direction cosines,momentum *
1253 /// * Output parameters *
1254 /// * VOUT Output co-ords,direction cosines,momentum *
1255 /// * User routine called *
1256 /// * CALL GUFLD(X,F) *
1257 /// * *
1258 /// * ==>Called by : <USER>, GUSWIM *
1259 /// * Authors R.Brun, M.Hansroul ********* *
1260 /// * V.Perevoztchikov (CUT STEP implementation) *
1261 /// * *
1262 /// * *
1263 /// ******************************************************************
1264
1265 Double_t h2, h4, f[4];
1266 Double_t /* xyzt[3], */ a, b, c, ph,ph2;
1267 Double_t secxs[4],secys[4],seczs[4],hxp[3];
1268 Double_t g1, g2, g3, g4, g5, g6, ang2, dxt, dyt, dzt;
1269 Double_t est, at, bt, ct, cba;
1270 Double_t f1, f2, f3, f4, rho, tet, hnorm, hp, rho1, sint, cost;
1271
1272 Double_t x;
1273 Double_t y;
1274 Double_t z;
1275
1276 Double_t xt;
1277 Double_t yt;
1278 Double_t zt;
1279
1280 // const Int_t maxit = 1992;
1281 const Int_t maxit = 500;
1282 const Int_t maxcut = 11;
1283
1284 const Double_t hmin = 1e-4; // !!! MT ADD, should be member
1285 const Double_t kdlt = 1e-3; // !!! MT CHANGE from 1e-4, should be member
1286 const Double_t kdlt32 = kdlt/32.;
1287 const Double_t kthird = 1./3.;
1288 const Double_t khalf = 0.5;
1289 const Double_t kec = 2.9979251e-3;
1290
1291 const Double_t kpisqua = 9.86960440109;
1292 const Int_t kix = 0;
1293 const Int_t kiy = 1;
1294 const Int_t kiz = 2;
1295 const Int_t kipx = 3;
1296 const Int_t kipy = 4;
1297 const Int_t kipz = 5;
1298
1299 // *.
1300 // *. ------------------------------------------------------------------
1301 // *.
1302 // * this constant is for units cm,gev/c and kgauss
1303 // *
1304 Int_t iter = 0;
1305 Int_t ncut = 0;
1306 for(Int_t j = 0; j < 7; j++)
1307 vout[j] = vect[j];
1308
1309 Double_t pinv = kec * fH.fCharge / vect[6];
1310 Double_t tl = 0.;
1311 Double_t h = step;
1312 Double_t rest;
1313
1314 do {
1315 rest = step - tl;
1316 if (TMath::Abs(h) > TMath::Abs(rest))
1317 h = rest;
1318
1319 f[0] = -fH.fB.fX;
1320 f[1] = -fH.fB.fY;
1321 f[2] = -fH.fB.fZ;
1322
1323 // * start of integration
1324 x = vout[0];
1325 y = vout[1];
1326 z = vout[2];
1327 a = vout[3];
1328 b = vout[4];
1329 c = vout[5];
1330
1331 h2 = khalf * h;
1332 h4 = khalf * h2;
1333 ph = pinv * h;
1334 ph2 = khalf * ph;
1335 secxs[0] = (b * f[2] - c * f[1]) * ph2;
1336 secys[0] = (c * f[0] - a * f[2]) * ph2;
1337 seczs[0] = (a * f[1] - b * f[0]) * ph2;
1338 ang2 = (secxs[0]*secxs[0] + secys[0]*secys[0] + seczs[0]*seczs[0]);
1339 if (ang2 > kpisqua) break;
1340
1341 dxt = h2 * a + h4 * secxs[0];
1342 dyt = h2 * b + h4 * secys[0];
1343 dzt = h2 * c + h4 * seczs[0];
1344 xt = x + dxt;
1345 yt = y + dyt;
1346 zt = z + dzt;
1347
1348 // * second intermediate point
1350 if (est > h) {
1351 if (ncut++ > maxcut) break;
1352 h *= khalf;
1353 continue;
1354 }
1355
1356 // xyzt[0] = xt;
1357 // xyzt[1] = yt;
1358 // xyzt[2] = zt;
1359
1361 f[0] = -fH.fB.fX;
1362 f[1] = -fH.fB.fY;
1363 f[2] = -fH.fB.fZ;
1364
1365 at = a + secxs[0];
1366 bt = b + secys[0];
1367 ct = c + seczs[0];
1368
1369 secxs[1] = (bt * f[2] - ct * f[1]) * ph2;
1370 secys[1] = (ct * f[0] - at * f[2]) * ph2;
1371 seczs[1] = (at * f[1] - bt * f[0]) * ph2;
1372 at = a + secxs[1];
1373 bt = b + secys[1];
1374 ct = c + seczs[1];
1375 secxs[2] = (bt * f[2] - ct * f[1]) * ph2;
1376 secys[2] = (ct * f[0] - at * f[2]) * ph2;
1377 seczs[2] = (at * f[1] - bt * f[0]) * ph2;
1378 dxt = h * (a + secxs[2]);
1379 dyt = h * (b + secys[2]);
1380 dzt = h * (c + seczs[2]);
1381 xt = x + dxt;
1382 yt = y + dyt;
1383 zt = z + dzt;
1384 at = a + 2.*secxs[2];
1385 bt = b + 2.*secys[2];
1386 ct = c + 2.*seczs[2];
1387
1389 if (est > 2.*TMath::Abs(h)) {
1390 if (ncut++ > maxcut) break;
1391 h *= khalf;
1392 continue;
1393 }
1394
1395 // xyzt[0] = xt;
1396 // xyzt[1] = yt;
1397 // xyzt[2] = zt;
1398
1400 f[0] = -fH.fB.fX;
1401 f[1] = -fH.fB.fY;
1402 f[2] = -fH.fB.fZ;
1403
1404 z = z + (c + (seczs[0] + seczs[1] + seczs[2]) * kthird) * h;
1405 y = y + (b + (secys[0] + secys[1] + secys[2]) * kthird) * h;
1406 x = x + (a + (secxs[0] + secxs[1] + secxs[2]) * kthird) * h;
1407
1408 secxs[3] = (bt*f[2] - ct*f[1])* ph2;
1409 secys[3] = (ct*f[0] - at*f[2])* ph2;
1410 seczs[3] = (at*f[1] - bt*f[0])* ph2;
1411 a = a+(secxs[0]+secxs[3]+2. * (secxs[1]+secxs[2])) * kthird;
1412 b = b+(secys[0]+secys[3]+2. * (secys[1]+secys[2])) * kthird;
1413 c = c+(seczs[0]+seczs[3]+2. * (seczs[1]+seczs[2])) * kthird;
1414
1415 est = TMath::Abs(secxs[0]+secxs[3] - (secxs[1]+secxs[2]))
1416 + TMath::Abs(secys[0]+secys[3] - (secys[1]+secys[2]))
1417 + TMath::Abs(seczs[0]+seczs[3] - (seczs[1]+seczs[2]));
1418
1419 if (est > kdlt && TMath::Abs(h) > hmin) {
1420 if (ncut++ > maxcut) break;
1421 h *= khalf;
1422 continue;
1423 }
1424
1425 ncut = 0;
1426 // * if too many iterations, go to helix
1427 if (iter++ > maxit) break;
1428
1429 tl += h;
1430 if (est < kdlt32)
1431 h *= 2.;
1432 cba = 1./ TMath::Sqrt(a*a + b*b + c*c);
1433 vout[0] = x;
1434 vout[1] = y;
1435 vout[2] = z;
1436 vout[3] = cba*a;
1437 vout[4] = cba*b;
1438 vout[5] = cba*c;
1439 rest = step - tl;
1440 if (step < 0.) rest = -rest;
1441 if (rest < 1.e-5*TMath::Abs(step))
1442 {
1443 Float_t dot = (vout[3]*vect[3] + vout[4]*vect[4] + vout[5]*vect[5]);
1444 fH.fPhi += TMath::ACos(dot);
1445 return;
1446 }
1447
1448 } while(true);
1449
1450 // angle too big, use helix
1451
1452 f1 = f[0];
1453 f2 = f[1];
1454 f3 = f[2];
1455 f4 = TMath::Sqrt(f1*f1+f2*f2+f3*f3);
1456 rho = -f4*pinv;
1457 tet = rho * step;
1458
1459 hnorm = 1./f4;
1460 f1 = f1*hnorm;
1461 f2 = f2*hnorm;
1462 f3 = f3*hnorm;
1463
1464 hxp[0] = f2*vect[kipz] - f3*vect[kipy];
1465 hxp[1] = f3*vect[kipx] - f1*vect[kipz];
1466 hxp[2] = f1*vect[kipy] - f2*vect[kipx];
1467
1468 hp = f1*vect[kipx] + f2*vect[kipy] + f3*vect[kipz];
1469
1470 rho1 = 1./rho;
1471 sint = TMath::Sin(tet);
1472 cost = 2.*TMath::Sin(khalf*tet)*TMath::Sin(khalf*tet);
1473
1474 g1 = sint*rho1;
1475 g2 = cost*rho1;
1476 g3 = (tet-sint) * hp*rho1;
1477 g4 = -cost;
1478 g5 = sint;
1479 g6 = cost * hp;
1480
1481 vout[kix] = vect[kix] + g1*vect[kipx] + g2*hxp[0] + g3*f1;
1482 vout[kiy] = vect[kiy] + g1*vect[kipy] + g2*hxp[1] + g3*f2;
1483 vout[kiz] = vect[kiz] + g1*vect[kipz] + g2*hxp[2] + g3*f3;
1484
1485 vout[kipx] = vect[kipx] + g4*vect[kipx] + g5*hxp[0] + g6*f1;
1486 vout[kipy] = vect[kipy] + g4*vect[kipy] + g5*hxp[1] + g6*f2;
1487 vout[kipz] = vect[kipz] + g4*vect[kipz] + g5*hxp[2] + g6*f3;
1488
1489 fH.fPhi += tet;
1490}
#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 a(i)
Definition RSha256.hxx:99
#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
unsigned char UChar_t
Unsigned Character 1 byte (unsigned char)
Definition RtypesCore.h:52
float Float_t
Float 4 bytes (float)
Definition RtypesCore.h:71
constexpr Bool_t kFALSE
Definition RtypesCore.h:108
double Double_t
Double 8 bytes.
Definition RtypesCore.h:73
constexpr Bool_t kTRUE
Definition RtypesCore.h:107
@ kRed
Definition Rtypes.h:67
@ kBlue
Definition Rtypes.h:67
@ kYellow
Definition Rtypes.h:67
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
Double_t Dot(const TGLVector3 &v1, const TGLVector3 &v2)
Definition TGLUtil.h:317
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t hmin
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t np
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
virtual void SetMarkerColor(Color_t mcolor=1)
Set the marker color.
Definition TAttMarker.h:39
virtual void SetMarkerStyle(Style_t mstyle=1)
Set the marker style.
Definition TAttMarker.h:41
virtual void SetMarkerSize(Size_t msize=1)
Set the marker size.
Definition TAttMarker.h:46
A list of TEveElements.
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.
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...
void Reset(Int_t n_points=0, Int_t n_int_ids=0)
Drop all data and set-up the data structures to recive new data.
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 ElementChanged(Bool_t update_scenes=kTRUE, Bool_t redraw=kFALSE) override
Element-change notification.
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
void LineToBounds(TEveVectorD &p)
Propagate neutral particle with momentum p to bounds.
Int_t GetCurrentPoint() const
Get index of current point on track.
static Double_t fgDefMagField
~TEveTrackPropagator() override
Destructor.
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.
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.
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.
void CheckReferenceCount(const TEveException &eh="TEveElement::CheckReferenceCount ") override
Check reference count - virtual from TEveElement.
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.
void OnZeroRefCount() override
Virtual from TEveRefBackPtr - track reference count has reached zero.
static const Double_t fgkB2C
void SetMaxStep(Double_t x)
Set maximum step-size and rebuild tracks.
Bool_t PointOverVertex(const TEveVector4D &v0, const TEveVector4D &v, Double_t *p=nullptr)
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.
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
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 Dot(const TEveVectorT &a) const
Definition TEveVector.h:168
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition TObject.cxx:1057
virtual Int_t SetNextPoint(Double_t x, Double_t y, Double_t z)
Set point following LastPoint to x, y, z.
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
TMath.
Definition TMathBase.h:35
Double_t ACos(Double_t)
Returns the principal value of the arc cosine of x, expressed in radians.
Definition TMath.h:643
Short_t Range(Short_t lb, Short_t ub, Short_t x)
Returns x if lb < x < up, lb if x < lb and ub if x > ub.
Definition TMathBase.h:303
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:673
Short_t Min(Short_t a, Short_t b)
Returns the smallest of a and b.
Definition TMathBase.h:199
Double_t Sin(Double_t)
Returns the sine of an angle of x radians.
Definition TMath.h:599
T * Cross(const T v1[3], const T v2[3], T out[3])
Calculates the Cross Product of two vectors: out = [v1 x v2].
Definition TMath.h:1284
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:124
constexpr Double_t TwoPi()
Definition TMath.h:47
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
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2339