From ccee46b404b908389e28d070d1593d376d38318b Mon Sep 17 00:00:00 2001
From: Matevz Tadel <matevz.tadel@cern.ch>
Date: Mon, 12 Jan 2009 19:04:45 +0000
Subject: [PATCH] From Alja.
EveTrackPropagator.cxx, TEveTrackPropagator.h:
Bug fix in track propagation in non-constant magnetic field.
Introduce maximum allowed step size.
tutorials/eve/track.C
Change macro to compare Helix and Runge Kutta stepper for
different types of magnetic field.
git-svn-id: http://root.cern.ch/svn/root/trunk@27119 27541ba8-7e3a-0410-8455-c3a389f83636
---
graf3d/eve/inc/TEveTrackPropagator.h | 28 ++--
graf3d/eve/src/TEveTrackPropagator.cxx | 189 ++++++++++++------------
tutorials/eve/track.C | 194 +++++++++++++++----------
3 files changed, 218 insertions(+), 193 deletions(-)
diff --git a/graf3d/eve/inc/TEveTrackPropagator.h b/graf3d/eve/inc/TEveTrackPropagator.h
index 626c9633bf7..e0a11bec467 100644
--- a/graf3d/eve/inc/TEveTrackPropagator.h
+++ b/graf3d/eve/inc/TEveTrackPropagator.h
@@ -45,8 +45,6 @@ public:
virtual TEveVector GetField(const TEveVector &v) const { return GetField(v.fX, v.fY, v.fZ);}
virtual TEveVector GetField(Float_t x, Float_t y, Float_t z) const = 0;
-
- virtual Float_t GetMaxFieldMag() const = 0;
ClassDef(TEveMagField, 0); // Abstract interface to magnetic field
};
@@ -69,8 +67,6 @@ public:
using TEveMagField::GetField;
virtual TEveVector GetField(Float_t /*x*/, Float_t /*y*/, Float_t /*z*/) const { return fB; }
- virtual Float_t GetMaxFieldMag() const { return fB.Mag(); }
-
ClassDef(TEveMagFieldConst, 0); // Interface to constant magnetic field.
};
@@ -87,7 +83,7 @@ protected:
Float_t fR2;
public:
- TEveMagFieldDuo(Float_t r, Float_t bIn, Float_t bOut) : TEveMagField(),
+ TEveMagFieldDuo(Float_t r, Float_t bIn, Float_t bOut) : TEveMagField(),
fBIn(0,0,bIn), fBOut(0,0,bOut), fR2(r*r)
{
fFieldConstant = kFALSE;
@@ -98,8 +94,6 @@ public:
virtual TEveVector GetField(Float_t x, Float_t y, Float_t /*z*/) const
{ return ((x*x+y*y)<fR2) ? fBIn : fBOut; }
- virtual Float_t GetMaxFieldMag() const { return TMath::Max(fBIn.Mag(), fBOut.Mag()); }
-
ClassDef(TEveMagFieldDuo, 0); // Interface to magnetic field with two different values depending of radius.
};
@@ -119,6 +113,8 @@ public:
Int_t fCharge; // Charge of tracked particle.
Float_t fMinAng; // Minimal angular step between two helix points.
Float_t fDelta; // Maximal error at the mid-point of the line connecting two helix points.
+ Float_t fMaxStep; // Maximum allowed step size.
+ Float_t fCurrentStep;
Float_t fPhi; // Accumulated angle to check fMaxOrbs by propagator.
Bool_t fValid; // Corner case pT~0 or B~0, possible in variable mag field.
@@ -137,18 +133,18 @@ public:
TEveVector fPt, fPl; // Transverse and longitudinal momentum.
Float_t fPtMag; // Magnitude of pT
Float_t fPlDir; // Momentum parallel to mag field.
- Float_t fTStep; // Transverse step arc-length in cm.
+ Float_t fLStep; // Transverse step arc-length in cm.
// ----------------------------------------------------------------
Helix_t();
- void UpdateHelix(const TEveVector & p, const TEveVector& b, Bool_t fullUpdate, Float_t fraction = -1);
- void UpdateRG (const TEveVector & p, const TEveVector& b, Float_t bMax = -1, Float_t maxStep = -1);
+ void UpdateHelix(const TEveVector & p, const TEveVector& b, Bool_t fullUpdate);
+ void UpdateRK (const TEveVector & p, const TEveVector& b);
void Step (const TEveVector4& v, const TEveVector& p, TEveVector4& vOut, TEveVector& pOut);
- Float_t GetStep() { return fTStep * TMath::Sqrt(1 + fLam*fLam); }
- Float_t GetStep2() { return fTStep * fTStep * (1 + fLam*fLam); }
+ Float_t GetStep() { return fLStep * TMath::Sqrt(1 + fLam*fLam); }
+ Float_t GetStep2() { return fLStep * fLStep * (1 + fLam*fLam); }
};
enum EStepper_e { kHelix, kRungeKutta };
@@ -157,7 +153,7 @@ private:
TEveTrackPropagator& operator=(const TEveTrackPropagator&); // Not implemented
protected:
- EStepper_e fStepper;
+ EStepper_e fStepper;
TEveMagField* fMagFieldObj;
@@ -168,8 +164,6 @@ protected:
// Helix limits
Float_t fMaxOrbs; // Maximal angular path of tracks' orbits (1 ~ 2Pi).
- Float_t fMaxStepRG; // Maximum step size in stepper RungeKuta.
-
// Path-mark / first-vertex control
Bool_t fEditPathMarks; // Show widgets for path-mark control in GUI editor.
Bool_t fFitDaughters; // Pass through daughter creation points when extrapolating a track.
@@ -192,7 +186,7 @@ protected:
Helix_t fH; // Helix.
void RebuildTracks();
- void Step(TEveVector4 &v, TEveVector &p, TEveVector4 &vOut, TEveVector &pOut, Float_t fraction);
+ void Step(TEveVector4 &v, TEveVector &p, TEveVector4 &vOut, TEveVector &pOut);
Bool_t LoopToVertex(TEveVector& v, TEveVector& p);
void LoopToBounds(TEveVector& p);
@@ -242,6 +236,7 @@ public:
void SetMaxOrbs(Float_t x);
void SetMinAng(Float_t x);
void SetDelta(Float_t x);
+ void SetMaxStep(Float_t x);
void SetEditPathMarks(Bool_t x) { fEditPathMarks = x; }
void SetRnrDaughters(Bool_t x);
@@ -264,6 +259,7 @@ public:
Float_t GetMaxOrbs() const { return fMaxOrbs; }
Float_t GetMinAng() const { return fH.fMinAng; }
Float_t GetDelta() const { return fH.fDelta; }
+ Float_t GetMaxStep() const { return fH.fMaxStep; }
Bool_t GetEditPathMarks() const { return fEditPathMarks; }
Bool_t GetRnrDaughters() const { return fRnrDaughters; }
diff --git a/graf3d/eve/src/TEveTrackPropagator.cxx b/graf3d/eve/src/TEveTrackPropagator.cxx
index 4c30211cf59..ff955637d34 100644
--- a/graf3d/eve/src/TEveTrackPropagator.cxx
+++ b/graf3d/eve/src/TEveTrackPropagator.cxx
@@ -16,19 +16,23 @@
#include <cassert>
+static Float_t kBMin = 1e-6;
+static Float_t kPtMinSqr = 1e-5;
+static Float_t kStepEps = 1e-3;
+
//______________________________________________________________________________
TEveTrackPropagator::Helix_t::Helix_t() :
fCharge(0), fMinAng(45), fDelta(0.1),
+ fMaxStep(20.f), fCurrentStep(20.f),
fPhi(0), fValid(kFALSE),
fLam(-1), fR(-1), fPhiStep(-1), fSin(-1), fCos(-1),
- fPtMag(-1), fPlDir(-1), fTStep(-1)
+ fPtMag(-1), fPlDir(-1), fLStep(-1)
{
// Default constructor.
}
//______________________________________________________________________________
-void TEveTrackPropagator::Helix_t::UpdateHelix(const TEveVector& p, const TEveVector& b,
- Bool_t fullUpdate, Float_t fraction)
+void TEveTrackPropagator::Helix_t::UpdateHelix(const TEveVector& p, const TEveVector& b, Bool_t fullUpdate)
{
// Update helix parameters.
@@ -52,34 +56,33 @@ void TEveTrackPropagator::Helix_t::UpdateHelix(const TEveVector& p, const TEveVe
if (fullUpdate)
{
using namespace TMath;
-
+ fValid = kFALSE;
Float_t a = fgkB2C * b.Mag() * Abs(fCharge);
- if (a > 1e-10 && fPtMag*fPtMag > 1e-12)
+ if (a > 1e-10 && fPtMag*fPtMag > kPtMinSqr)
{
fValid = kTRUE;
fR = Abs(fPtMag / a);
fLam = fPl.Mag() / fPtMag;
if (fPlDir < 0) fLam = - fLam;
- }
- else
- {
- fValid = kFALSE;
- return;
- }
- // phi steps
- fPhiStep = fMinAng * DegToRad();
- if (fDelta < fR)
- {
- Float_t ang = 2*ACos(1 - fDelta/fR);
- if (ang < fPhiStep) fPhiStep = ang;
- }
- if (fraction > 0) fPhiStep *= fraction;
+ // get phi step, compare fMinAng with fDelta
+ fPhiStep = fMinAng * DegToRad();
+ if (fDelta < fR)
+ {
+ Float_t ang = 2*ACos(1 - fDelta/fR);
+ if (ang < fPhiStep) fPhiStep = ang;
+ }
+
+ // check max step size
+ fCurrentStep = fR*fPhiStep*(1 + fLam*fLam);
+ if (fCurrentStep > fMaxStep)
+ fPhiStep *= (fMaxStep/fCurrentStep);
- fTStep = fR*fPhiStep;
- fSin = Sin(fPhiStep);
- fCos = Cos(fPhiStep);
+ fLStep = fR*fPhiStep*fLam;
+ fSin = Sin(fPhiStep);
+ fCos = Cos(fPhiStep);
+ }
}
}
@@ -93,9 +96,9 @@ void TEveTrackPropagator::Helix_t::Step(const TEveVector4& v, const TEveVector&
if (fValid)
{
- TEveVector d = fE2*(fR*fSin) + fE3*(fR*(1-fCos)) + fE1*(fLam*fTStep);
+ TEveVector d = fE2*(fR*fSin) + fE3*(fR*(1-fCos)) + fE1*fLStep;
vOut += d;
- vOut.fT += fTStep;
+ vOut.fT += fLStep;
pOut = fPl + fE2*(fPtMag*fCos) + fE3*(fPtMag*fSin);
@@ -103,44 +106,31 @@ void TEveTrackPropagator::Helix_t::Step(const TEveVector4& v, const TEveVector&
}
else
{
- // case: pT < 1e-6 or B < 1e-7
+ // case: pT < kPtMinSqr or B < kBMin
// might happen if field directon changes pT ~ 0 or B becomes zero
-
- if (fTStep == -1)
- {
- printf("WARNING TEveTrackPropagator::Helix_t::Step step-size not initialised.\n");
- }
- vOut += p * (fTStep / p.Mag());
+ vOut += p * (fMaxStep / p.Mag());
pOut = p;
}
}
//______________________________________________________________________________
-void TEveTrackPropagator::Helix_t::UpdateRG(const TEveVector& p, const TEveVector& b, Float_t bMax, Float_t maxStep)
+void TEveTrackPropagator::Helix_t::UpdateRK(const TEveVector& p, const TEveVector& b)
{
- // Update helix for stepper RungeKutta.
-
- if (bMax > 0)
- {
- // full update
- Float_t a = fgkB2C * bMax * fCharge;
- fPhiStep = fMinAng * TMath::DegToRad();
- if (a)
- {
- fR = TMath::Abs(p.Mag() / a);
- fTStep = TMath::Min(fR*fPhiStep, maxStep);
+ // Update helix for stepper RungeKutta.
+
+ if (TMath::Abs(b.Mag()) > kBMin && fCharge)
+ {
fValid = kTRUE;
- }
- else
- {
- fValid = kFALSE;
- fTStep = maxStep;
- }
- }
- fB = b;
- fPlDir = p.Dot(fB);
+ // cached values for propagator
+ fB = b;
+ fPlDir = p.Dot(fB);
+ }
+ else
+ {
+ fValid = kFALSE;
+ }
}
@@ -179,7 +169,6 @@ TEveTrackPropagator::TEveTrackPropagator(const Text_t* n, const Text_t* t,
fNMax (4096),
fMaxOrbs (0.5),
- fMaxStepRG (10),
fEditPathMarks (kTRUE),
fFitDaughters (kTRUE),
@@ -295,10 +284,9 @@ Bool_t TEveTrackPropagator::GoToVertex(TEveVector& v, TEveVector& p)
if (fStepper == kHelix)
fH.UpdateHelix(p, fMagFieldObj->GetField(fV), kTRUE);
else
- fH.UpdateRG(p, fMagFieldObj->GetField(fV), fMagFieldObj->GetMaxFieldMag(), fMaxStepRG);
+ fH.UpdateRK(p, fMagFieldObj->GetField(fV));
-
- if ((v-fV).Mag() < 1e-3)
+ if ((v-fV).Mag() < kStepEps)
{
fPoints.push_back(v);
return kTRUE;
@@ -312,55 +300,50 @@ void TEveTrackPropagator::GoToBounds(TEveVector& p)
{
// Propagate particle to bounds.
// Return TRUE if hit bounds.
-
+
if (fStepper == kHelix)
fH.UpdateHelix(p, fMagFieldObj->GetField(fV), kTRUE);
else
- fH.UpdateRG(p, fMagFieldObj->GetField(fV), fMagFieldObj->GetMaxFieldMag(), fMaxStepRG);
+ fH.UpdateRK(p, fMagFieldObj->GetField(fV));
-
fH.fValid ? LoopToBounds(p): LineToBounds(p);
}
//______________________________________________________________________________
-void TEveTrackPropagator::Step(TEveVector4 &v, TEveVector &p, TEveVector4 &vOut, TEveVector &pOut,
- Float_t fraction = -1)
+void TEveTrackPropagator::Step(TEveVector4 &v, TEveVector &p, TEveVector4 &vOut, TEveVector &pOut)
{
// Wrapper to step helix.
if (fStepper == kHelix)
{
- fH.UpdateHelix(p, fH.fB, !fMagFieldObj->IsConst(), fraction);
+ fH.UpdateHelix(p, fMagFieldObj->GetField(v), !fMagFieldObj->IsConst());
fH.Step(v, p, vOut, pOut);
}
else
{
- fH.UpdateRG(p, fMagFieldObj->GetField(fV));
- Float_t step = fH.fTStep;
- if (fraction > 0)
- step *= fraction;
-
- Double_t vecRGI[7];
- vecRGI[0] = v.fX;
- vecRGI[1] = v.fY;
- vecRGI[2] = v.fZ;
+ fH.UpdateRK(p, fMagFieldObj->GetField(v));
+
+ Double_t vecRKIn[7];
+ vecRKIn[0] = v.fX;
+ vecRKIn[1] = v.fY;
+ vecRKIn[2] = v.fZ;
Float_t nm = 1/p.Mag();
- vecRGI[3] = p.fX*nm;
- vecRGI[4] = p.fY*nm;
- vecRGI[5] = p.fZ*nm;
- vecRGI[6] = p.Mag();
-
- Double_t vecRGO[7];
- OneStepRungeKutta(fH.fCharge, step, vecRGI, vecRGO);
-
- vOut.fX = vecRGO[0];
- vOut.fY = vecRGO[1];
- vOut.fZ = vecRGO[2];
- vOut.fT += step;
- Double_t pm = vecRGO[6];
- pOut.fX = vecRGO[3]*pm;
- pOut.fY = vecRGO[4]*pm;
- pOut.fZ = vecRGO[5]*pm;
+ vecRKIn[3] = p.fX*nm;
+ vecRKIn[4] = p.fY*nm;
+ vecRKIn[5] = p.fZ*nm;
+ vecRKIn[6] = p.Mag();
+
+ Double_t vecRKOut[7];
+ OneStepRungeKutta(fH.fCharge, fH.fMaxStep, vecRKIn, vecRKOut);
+
+ vOut.fX = vecRKOut[0];
+ vOut.fY = vecRKOut[1];
+ vOut.fZ = vecRKOut[2];
+ vOut.fT += fH.fMaxStep;
+ Double_t pm = vecRKOut[6];
+ pOut.fX = vecRKOut[3]*pm;
+ pOut.fY = vecRKOut[4]*pm;
+ pOut.fZ = vecRKOut[5]*pm;
}
}
@@ -468,22 +451,23 @@ Bool_t TEveTrackPropagator::LoopToVertex(TEveVector& v, TEveVector& p)
TEveVector pln = fH.fPl;
pln.Normalize();
- if (pln.Dot(d2) > 1e-10)
+ if (d1.Mag() > kStepEps)
{
- // set step fraction
- Float_t frac = pln.Dot(d1)/pln.Dot(d2);
- // printf("frac %f %d p (%f %f %f)\n", frac, np-first_point, fH.fPl.fX, fH.fPl.fY, fH.fPl.fZ);
-
+ // step for forced step size;
+ Float_t origMaxStep = fH.fMaxStep;
+ fH.fMaxStep = d1.Mag();
if (fStepper == kHelix)
- fH.UpdateHelix(p, fH.fB, kTRUE, frac);
+ fH.UpdateHelix(p, fMagFieldObj->GetField(fV), kTRUE);
else
- fH.fTStep *= frac;
+ fH.UpdateRK(p, fMagFieldObj->GetField(fV));
Step(currV, p, forwV, forwP);
p = forwP;
currV = forwV;
fPoints.push_back(currV); np++;
+ fH.fMaxStep = origMaxStep;
+ // distibute offset
TEveVector off(v); off -= currV;
off *= 1.0f / currV.fT;
for (Int_t i = first_point; i < np; ++i)
@@ -561,7 +545,7 @@ Bool_t TEveTrackPropagator::HelixIntersectPlane(const TEveVector& p,
TEveVector pos(fV);
TEveVector mom(p);
if (fMagFieldObj->IsConst())
- fH.UpdateHelix(mom, fMagFieldObj->GetField(pos), fH.fCharge, kTRUE);
+ fH.UpdateHelix(mom, fMagFieldObj->GetField(pos), kFALSE);
TEveVector n(normal);
TEveVector delta = pos - point;
@@ -636,7 +620,7 @@ Bool_t TEveTrackPropagator::IntersectPlane(const TEveVector& p,
// Returns:
// kFALSE if intersection can not be found, kTRUE otherwise.
- if (fH.fCharge && fMagFieldObj && p.Perp2() > 1e-12)
+ if (fH.fCharge && fMagFieldObj && p.Perp2() > kPtMinSqr)
return HelixIntersectPlane(p, point, normal, itsect);
else
return LineIntersectPlane(p, point, normal, itsect);
@@ -818,6 +802,14 @@ void TEveTrackPropagator::SetRnrReferences(Bool_t rnr)
RebuildTracks();
}
+//______________________________________________________________________________
+void TEveTrackPropagator::SetMaxStep(Float_t x)
+{
+ // Set maximum radius and rebuild tracks.
+
+ fH.fMaxStep = x;
+ RebuildTracks();
+}
//______________________________________________________________________________
void TEveTrackPropagator::OneStepRungeKutta(Double_t charge, Double_t step,
@@ -897,7 +889,6 @@ void TEveTrackPropagator::OneStepRungeKutta(Double_t charge, Double_t step,
Double_t h = step;
Double_t rest;
-
do {
rest = step - tl;
if (TMath::Abs(h) > TMath::Abs(rest)) h = rest;
@@ -950,6 +941,7 @@ void TEveTrackPropagator::OneStepRungeKutta(Double_t charge, Double_t step,
xyzt[2] = zt;
//cmodif: call gufld(xyzt,f) changed into:
+ fH.fB = fMagFieldObj->GetField(xt, yt, zt);
f[0] = -fH.fB.fX;
f[1] = -fH.fB.fY;
f[2] = -fH.fB.fZ;
@@ -989,6 +981,7 @@ void TEveTrackPropagator::OneStepRungeKutta(Double_t charge, Double_t step,
xyzt[2] = zt;
//cmodif: call gufld(xyzt,f) changed into:
+ fH.fB = fMagFieldObj->GetField(xt, yt, zt);
f[0] = -fH.fB.fX;
f[1] = -fH.fB.fY;
f[2] = -fH.fB.fZ;
@@ -1074,7 +1067,7 @@ void TEveTrackPropagator::OneStepRungeKutta(Double_t charge, Double_t step,
vout[kipx] = vect[kipx] + g4*vect[kipx] + g5*hxp[0] + g6*f1;
vout[kipy] = vect[kipy] + g4*vect[kipy] + g5*hxp[1] + g6*f2;
- vout[kipz] = vect[kipz] + g4*vect[kipz] + g5*hxp[2] + g6*f3;
+ vout[kipz] = vect[kipz] + g4*vect[kipz] + g5*hxp[2] + g6*f3;
fH.fPhi += tet;
}
diff --git a/tutorials/eve/track.C b/tutorials/eve/track.C
index 7859e86fb4a..251e798a5e8 100644
--- a/tutorials/eve/track.C
+++ b/tutorials/eve/track.C
@@ -3,55 +3,31 @@
// Makes some tracks with three different magnetic field types.
-void track(Int_t bCase = 2)
+#include "TEveTrackPropagator.h"
+#include "TEveTrack.h"
+#include "TEveManager.h"
+#include "TEveViewer.h"
+#include "TSystem.h"
+#include "TGLViewer.h"
+#include "TMath.h"
+
+#include "TEveViewer.h"
+#include "TEvePointSet.h"
+
+class GappedField : public TEveMagField
{
- gSystem->IgnoreSignal(kSigSegmentationViolation, true);
- TEveManager::Create();
-
- TEveTrackPropagator* prop = new TEveTrackPropagator();
- prop->SetRnrDaughters(kTRUE);
- gEve->AddElement(prop);
-
- TEveTrack *track = 0;
- switch (bCase)
- {
- case 0:
- {
- // B = 0 no difference btween signed and charge particles
- prop->SetMagField(0.);
- prop->SetName("ZeroB");
- track = make_track(prop, 0);
- track = make_track(prop, 1);
- break;
- }
-
- case 1:
- {
- // constant B field, const angle
- prop->SetMagFieldObj(new TEveMagFieldConst(0., 0., -3.8));
- prop->SetName("ConstB");
- prop->SetMaxR(123);
- prop->SetMaxZ(300);
- track = make_track(prop, 0);
- track = make_track(prop, -1);
- break;
- }
- case 2:
- {
- // variable B field, sign change at R = 200 cm
- prop->SetMagFieldObj(new TEveMagFieldDuo(200, -4.4, 2));
- prop->SetName("DuoB");
- prop->SetMaxR(600);
- prop->SetMaxZ(700);
- prop->SetMinAng(1.f);
- track = make_track(prop, 0);
- track = make_track(prop, -1);
- break;
- }
- };
- // make_projected(track);
- gEve->Redraw3D(1);
-}
+public:
+ GappedField():TEveMagField(){}
+ ~GappedField(){};
+ using TEveMagField::GetField;
+
+ virtual TEveVector GetField(Float_t /*x*/, Float_t /*y*/, Float_t z) const
+ {
+ if (TMath::Abs(z) < 300) return TEveVector(0, 0, -4);
+ if (TMath::Abs(z) < 600) return TEveVector(0, 0, 0);
+ return TEveVector(0, 0, 4);
+ }
+};
//______________________________________________________________________________
TEveTrack* make_track(TEveTrackPropagator* prop, Int_t sign)
@@ -63,9 +39,10 @@ TEveTrack* make_track(TEveTrackPropagator* prop, Int_t sign)
rc->fV.Set(0.028558, -0.000918, 3.691919);
rc->fP.Set(0.767095, -2.400006, -0.313103);
rc->fSign = sign;
+
+
TEveTrack* track = new TEveTrack(rc, prop);
track->SetName(Form("Charge %d", sign));
-
// daughter 0
TEvePathMark* pm1 = new TEvePathMark(TEvePathMark::kDaughter);
pm1->fV.Set(1.479084, -4.370661, 3.119761);
@@ -75,38 +52,97 @@ TEveTrack* make_track(TEveTrackPropagator* prop, Int_t sign)
pm2->fV.Set(57.72345, -89.77011, -9.783746);
track->AddPathMark(*pm2);
- // color by sign
- if (sign)
- {
- track->SetLineColor(kRed-3);
- track->SetLineWidth(2);
- }
- else
- {
- track->SetLineColor(kOrange);
- }
-
- track->MakeTrack();
- gEve->AddElement(track, prop);
-
return track;
}
-//______________________________________________________________________________
-void make_projected(TEveTrack *track)
+
+void track(Int_t bCase = 3, Bool_t isRungeKutta = kTRUE)
{
- // Make non-linear projected view.
-
- TEveViewer* v1 = gEve->SpawnNewViewer("2D Viewer");
- TEveScene* s1 = gEve->SpawnNewScene("Projected Event");
- v1->AddScene(s1);
- TGLViewer* v = v1->GetGLViewer();
- v->SetCurrentCamera(TGLViewer::kCameraOrthoXOY);
- v->SetGuideState(TGLUtil::kAxesOrigin, kTRUE, kFALSE, 0);
-
- TEveProjectionManager* mng = new TEveProjectionManager();
- mng->SetProjection(TEveProjection::kPT_RhoZ);
- gEve->AddElement(mng, s1);
- gEve->AddToListTree(mng, kTRUE);
- mng->ImportElements(track);
+#if defined (__CINT__)
+ Error("track.C", "Must be run in compiled mode!");
+ return;
+#endif
+
+ gSystem->IgnoreSignal(kSigSegmentationViolation, true);
+ TEveManager::Create();
+
+ TEveTrackList *list = new TEveTrackList();
+ TEveTrackPropagator* prop = list->GetPropagator();
+ prop->SetFitDaughters(kFALSE);
+ prop->SetMaxZ(1000);
+
+ if (isRungeKutta)
+ {
+ prop->SetStepper(TEveTrackPropagator::kRungeKutta);
+ list->SetName("RK Propagator");
+ }
+ else
+ {
+ list->SetName("Heix Propagator");
+ }
+
+ TEveTrack *track = 0;
+ switch (bCase)
+ {
+ case 0:
+ {
+ // B = 0 no difference btween signed and charge particles
+ prop->SetMagField(0.);
+ list->SetElementName(Form("%s, zeroB", list->GetElementName()));
+ track = make_track(prop, 1);
+ break;
+ }
+
+ case 1:
+ {
+ // constant B field, const angle
+ prop->SetMagFieldObj(new TEveMagFieldConst(0., 0., -3.8));
+ list->SetElementName(Form("%s, constB", list->GetElementName()));
+ track = make_track(prop, 1);
+ break;
+ }
+ case 2:
+ {
+ // variable B field, sign change at R = 200 cm
+ prop->SetMagFieldObj(new TEveMagFieldDuo(200, -4.4, 2));
+ list->SetElementName(Form("%s, duoB", list->GetElementName()));
+ track = make_track(prop, 1);
+ break;
+ }
+ case 3:
+ {
+ // gapped field
+ prop->SetMagFieldObj(new GappedField());
+ list->SetElementName(Form("%s, gappedB", list->GetElementName()));
+
+
+ TEveRecTrack *rc = new TEveRecTrack();
+ rc->fV.Set(0.028558, -0.000918, 3.691919);
+ rc->fP.Set(0.767095, -0.400006, 2.313103);
+ rc->fSign = 1;
+ track = new TEveTrack(rc, prop);
+
+ TEvePointSet* marker = new TEvePointSet(2);
+ marker->SetElementName("B field break points");
+ marker->SetPoint(0, 0., 0., 300.f);
+ marker->SetPoint(1, 0., 0., 600.f);
+ marker->SetMarkerColor(3);
+ gEve->AddElement(marker);
+ }
+ };
+
+ if (isRungeKutta)
+ list->SetLineColor(kMagenta);
+ else
+ list->SetLineColor(kCyan);
+
+ track->SetLineColor(list->GetLineColor());
+
+ gEve->AddElement(track, list);
+ gEve->AddElement(list);
+ track->MakeTrack();
+
+ TEveViewer* v = gEve->GetDefaultViewer();
+ v->GetGLViewer()->SetGuideState(TGLUtil::kAxesOrigin, kTRUE, kFALSE, 0);
+ gEve->Redraw3D(1);
}
--
GitLab