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