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