Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
track.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_eve
3 /// Demonstrates usage of TEveTrackPRopagator with different magnetic
4 /// field configurations.
5 /// Needs to be run in compiled mode.
6 /// root
7 /// ~~~{.cpp}
8 /// .L track.C+
9 /// track(3, kTRUE)
10 /// ~~~
11 /// void track(Int_t mode = 5, Bool_t isRungeKutta = kTRUE)
12 /// Modes are
13 /// 0. B = 0, no difference between signed and charge particles;
14 /// 1. constant B field (along z, but could have arbitrary direction);
15 /// 2. variable B field, sign change at R = 200 cm;
16 /// 3. magnetic field with a zero-field region;
17 /// 4. CMS magnetic field - simple track;
18 /// 5. CMS magnetic field - track with different path-marks.
19 /// 6. Conceptual ILC detector, problematic track
20 ///
21 /// \image html eve_track.png
22 /// \macro_code
23 ///
24 /// \author Alja Mrak-Tadel
25 
26 
27 #include "TEveTrackPropagator.h"
28 #include "TEveTrack.h"
29 #include "TEveVSDStructs.h"
30 #include "TEveManager.h"
31 #include "TEveViewer.h"
32 #include "TSystem.h"
33 #include "TGLViewer.h"
34 #include "TMath.h"
35 
36 #include "TEveViewer.h"
37 #include "TEvePointSet.h"
38 
39 #include <iostream>
40 
41 TEveTrackPropagator* g_prop = 0;
42 
43 class GappedField : public TEveMagField
44 {
45 public:
46  GappedField():TEveMagField(){}
47  ~GappedField(){};
48  using TEveMagField::GetField;
49 
50  virtual TEveVectorD GetFieldD(Double_t /*x*/, Double_t /*y*/, Double_t z) const
51  {
52  if (TMath::Abs(z) < 300) return TEveVectorD(0, 0, -4);
53  if (TMath::Abs(z) < 600) return TEveVectorD(0, 0, 0);
54  return TEveVectorD(0, 0, 4);
55  }
56 };
57 
58 //==============================================================================
59 
60 class CmsMagField: public TEveMagField
61 {
62  bool m_magnetIsOn;
63  bool m_reverse;
64  bool m_simpleModel;
65 
66 public:
67  CmsMagField():
68  m_magnetIsOn(true),
69  m_reverse(false),
70  m_simpleModel(true){}
71 
72  virtual ~CmsMagField(){}
73  virtual Double_t GetMaxFieldMagD() const { return m_magnetIsOn ? 3.8 : 0.0; }
74  void setMagnetState( bool state )
75  {
76  if (state != m_magnetIsOn)
77  {
78  if ( state )
79  std::cout << "Magnet state is changed to ON" << std::endl;
80  else
81  std::cout << "Magnet state is changed to OFF" << std::endl;
82  }
83  m_magnetIsOn = state;
84  }
85 
86  bool isMagnetOn() const { return m_magnetIsOn;}
87  void setReverseState(bool state) { m_reverse = state; }
88  bool isReverse() const { return m_reverse;}
89  void setSimpleModel(bool simpleModel) { m_simpleModel = simpleModel; }
90  bool isSimpleModel() const { return m_simpleModel;}
91 
92  using TEveMagField::GetField;
93 
94  virtual TEveVectorD GetFieldD(Double_t x, Double_t y, Double_t z) const
95  {
96  double R = sqrt(x*x+y*y);
97  double field = m_reverse?-GetMaxFieldMag():GetMaxFieldMag();
98  //barrel
99  if ( TMath::Abs(z)<724 )
100  {
101  //inside solenoid
102  if ( R < 300) return TEveVectorD(0,0,field);
103  // outside solenoid
104  if ( m_simpleModel ||
105  ( R>461.0 && R<490.5 ) ||
106  ( R>534.5 && R<597.5 ) ||
107  ( R>637.0 && R<700.0 ) )
108  return TEveVectorD(0,0,-field/3.8*1.2);
109 
110  } else {
111  // endcaps
112  if (m_simpleModel)
113  {
114  if ( R < 50 ) return TEveVectorD(0,0,field);
115  if ( z > 0 )
116  return TEveVectorD(x/R*field/3.8*2.0, y/R*field/3.8*2.0, 0);
117  else
118  return TEveVectorD(-x/R*field/3.8*2.0, -y/R*field/3.8*2.0, 0);
119  }
120  // proper model
121  if ( ( TMath::Abs(z)>724 && TMath::Abs(z)<786 ) ||
122  ( TMath::Abs(z)>850 && TMath::Abs(z)<910 ) ||
123  ( TMath::Abs(z)>975 && TMath::Abs(z)<1003 ) )
124  {
125  if ( z > 0 )
126  return TEveVectorD(x/R*field/3.8*2.0, y/R*field/3.8*2.0, 0);
127  else
128  return TEveVectorD(-x/R*field/3.8*2.0, -y/R*field/3.8*2.0, 0);
129  }
130  }
131  return TEveVectorD(0,0,0);
132  }
133 };
134 
135 
136 //==============================================================================
137 //==============================================================================
138 
139 //______________________________________________________________________________
140 TEveTrack* make_track(TEveTrackPropagator* prop, Int_t sign)
141 {
142  // Make track with given propagator.
143  // Add to math-marks to test fit.
144 
145  auto rc = new TEveRecTrackD();
146  rc->fV.Set(0.028558, -0.000918, 3.691919);
147  rc->fP.Set(0.767095, -2.400006, -0.313103);
148  rc->fSign = sign;
149 
150  auto track = new TEveTrack(rc, prop);
151  track->SetName(Form("Charge %d", sign));
152  // daughter 0
153  auto pm1 = new TEvePathMarkD(TEvePathMarkD::kDaughter);
154  pm1->fV.Set(1.479084, -4.370661, 3.119761);
155  track->AddPathMark(*pm1);
156  // daughter 1
157  auto pm2 = new TEvePathMarkD(TEvePathMarkD::kDaughter);
158  pm2->fV.Set(57.72345, -89.77011, -9.783746);
159  track->AddPathMark(*pm2);
160 
161  return track;
162 }
163 
164 
165 void track(Int_t mode = 1, Bool_t isRungeKutta = kTRUE)
166 {
167 
168  gSystem->IgnoreSignal(kSigSegmentationViolation, true);
169  TEveManager::Create();
170 
171  auto list = new TEveTrackList();
172  auto prop = g_prop = list->GetPropagator();
173  prop->SetFitDaughters(kFALSE);
174  prop->SetMaxZ(1000);
175 
176  if (isRungeKutta)
177  {
178  prop->SetStepper(TEveTrackPropagator::kRungeKutta);
179  list->SetName("RK Propagator");
180  }
181  else
182  {
183  list->SetName("Heix Propagator");
184  }
185 
186  TEveTrack *track = 0;
187  switch (mode)
188  {
189  case 0:
190  {
191  // B = 0 no difference between signed and charge particles
192  prop->SetMagField(0.);
193  list->SetElementName(Form("%s, zeroB", list->GetElementName()));
194  track = make_track(prop, 1);
195  break;
196  }
197 
198  case 1:
199  {
200  // constant B field, const angle
201  prop->SetMagFieldObj(new TEveMagFieldConst(0., 0., -3.8));
202  list->SetElementName(Form("%s, constB", list->GetElementName()));
203  track = make_track(prop, 1);
204  break;
205  }
206 
207  case 2:
208  {
209  // variable B field, sign change at R = 200 cm
210  prop->SetMagFieldObj(new TEveMagFieldDuo(200, -4.4, 2));
211  list->SetElementName(Form("%s, duoB", list->GetElementName()));
212  track = make_track(prop, 1);
213  break;
214  }
215 
216  case 3:
217  {
218  // gapped field
219  prop->SetMagFieldObj(new GappedField());
220  list->SetElementName(Form("%s, gappedB", list->GetElementName()));
221 
222 
223  auto rc = new TEveRecTrackD();
224  rc->fV.Set(0.028558, -0.000918, 3.691919);
225  rc->fP.Set(0.767095, -0.400006, 2.313103);
226  rc->fSign = 1;
227  track = new TEveTrack(rc, prop);
228 
229  auto marker = new TEvePointSet(2);
230  marker->SetElementName("B field break points");
231  marker->SetPoint(0, 0., 0., 300.f);
232  marker->SetPoint(1, 0., 0., 600.f);
233  marker->SetMarkerColor(3);
234  gEve->AddElement(marker);
235  break;
236  }
237 
238  case 4:
239  {
240  // Magnetic field of CMS I.
241  auto mf = new CmsMagField;
242  mf->setReverseState(true);
243 
244  prop->SetMagFieldObj(mf);
245  prop->SetMaxR(1000);
246  prop->SetMaxZ(1000);
247  prop->SetRnrDaughters(kTRUE);
248  prop->SetRnrDecay(kTRUE);
249  prop->RefPMAtt().SetMarkerStyle(4);
250  list->SetElementName(Form("%s, CMS field", list->GetElementName()));
251 
252 
253  auto rc = new TEveRecTrackD();
254  rc->fV.Set(0.027667, 0.007919, 0.895964);
255  rc->fP.Set(3.903134, 2.252232, -3.731366);
256  rc->fSign = -1;
257  track = new TEveTrack(rc, prop);
258 
259  track->AddPathMark(TEvePathMarkD(TEvePathMarkD::kDaughter,
260  TEveVectorD(3.576755e+00, 2.080579e+00, -2.507230e+00)));
261  track->AddPathMark(TEvePathMarkD(TEvePathMarkD::kDaughter,
262  TEveVectorD(8.440379e+01, 6.548286e+01, -8.788129e+01)));
263  track->AddPathMark(TEvePathMarkD(TEvePathMarkD::kDaughter,
264  TEveVectorD(1.841321e+02, 3.915693e+02, -3.843072e+02)));
265  track->AddPathMark(TEvePathMarkD(TEvePathMarkD::kDaughter,
266  TEveVectorD(1.946167e+02, 4.793932e+02, -4.615060e+02)));
267  track->AddPathMark(TEvePathMarkD(TEvePathMarkD::kDecay,
268  TEveVectorD(2.249656e+02, 5.835767e+02, -5.565275e+02)));
269 
270  track->SetRnrPoints(kTRUE);
271  track->SetMarkerStyle(4);
272 
273  break;
274  }
275 
276  case 5:
277  {
278  // Magnetic field of CMS I.
279  auto mf = new CmsMagField;
280  mf->setReverseState(true);
281  mf->setSimpleModel(false);
282 
283  prop->SetMagFieldObj(mf);
284  prop->SetMaxR(1000);
285  prop->SetMaxZ(1000);
286  prop->SetRnrReferences(kTRUE);
287  prop->SetRnrDaughters(kTRUE);
288  prop->SetRnrDecay(kTRUE);
289  prop->RefPMAtt().SetMarkerStyle(4);
290  list->SetElementName(Form("%s, CMS field", list->GetElementName()));
291 
292  auto rc = new TEveRecTrackD();
293  rc->fV.Set(-16.426592, 16.403185, -19.782692);
294  rc->fP.Set(3.631100, 3.643450, 0.682254);
295  rc->fSign = -1;
296  track = new TEveTrack(rc, prop);
297 
298  track->AddPathMark(TEvePathMarkD(TEvePathMarkD::kReference,
299  TEveVectorD(-1.642659e+01, 1.640318e+01, -1.978269e+01),
300  TEveVectorD(3.631100, 3.643450, 0.682254)));
301  track->AddPathMark(TEvePathMarkD(TEvePathMarkD::kReference,
302  TEveVectorD(-1.859987e+00, 3.172243e+01, -1.697866e+01),
303  TEveVectorD(3.456056, 3.809894, 0.682254)));
304  track->AddPathMark(TEvePathMarkD(TEvePathMarkD::kReference,
305  TEveVectorD(4.847579e+01, 9.871711e+01, -5.835719e+00),
306  TEveVectorD(2.711614, 4.409945, 0.687656)));
307  track->AddPathMark(TEvePathMarkD(TEvePathMarkD::kDaughter,
308  TEveVectorD(1.342045e+02, 4.203950e+02, 3.846268e+01)));
309  track->AddPathMark(TEvePathMarkD(TEvePathMarkD::kDaughter,
310  TEveVectorD(1.483827e+02, 5.124750e+02, 5.064311e+01)));
311  track->AddPathMark(TEvePathMarkD(TEvePathMarkD::kDaughter,
312  TEveVectorD(1.674676e+02, 6.167731e+02, 6.517403e+01)));
313  track->AddPathMark(TEvePathMarkD(TEvePathMarkD::kDecay,
314  TEveVectorD(1.884976e+02, 7.202000e+02, 7.919290e+01)));
315 
316  track->SetRnrPoints(kTRUE);
317  track->SetMarkerStyle(4);
318 
319  break;
320  }
321 
322  case 6:
323  {
324  // Problematic track from Druid
325  prop->SetMagFieldObj(new TEveMagFieldDuo(350, -3.5, 2.0));
326  prop->SetMaxR(1000);
327  prop->SetMaxZ(1000);
328  prop->SetRnrReferences(kTRUE);
329  prop->SetRnrDaughters(kTRUE);
330  prop->SetRnrDecay(kTRUE);
331  prop->RefPMAtt().SetMarkerStyle(4);
332  list->SetElementName(Form("%s, Some ILC Detector field",
333  list->GetElementName()));
334 
335  auto rc = new TEveRecTrackD();
336  rc->fV.Set(57.1068, 31.2401, -7.07629);
337  rc->fP.Set(4.82895, 2.35083, -0.611757);
338  rc->fSign = 1;
339  track = new TEveTrack(rc, prop);
340 
341  track->AddPathMark(TEvePathMarkD(TEvePathMarkD::kDaughter,
342  TEveVectorD(1.692235e+02, 7.047929e+01, -2.064785e+01)));
343  track->AddPathMark(TEvePathMarkD(TEvePathMarkD::kDaughter,
344  TEveVectorD(5.806180e+02, 6.990633e+01, -6.450000e+01)));
345  track->AddPathMark(TEvePathMarkD(TEvePathMarkD::kDecay,
346  TEveVectorD(6.527213e+02, 1.473249e+02, -8.348498e+01)));
347 
348  track->SetRnrPoints(kTRUE);
349  track->SetMarkerStyle(4);
350 
351  break;
352  }
353  };
354 
355  if (isRungeKutta)
356  list->SetLineColor(kMagenta);
357  else
358  list->SetLineColor(kCyan);
359 
360  track->SetLineColor(list->GetLineColor());
361 
362  gEve->AddElement(list);
363  list->AddElement(track);
364 
365  track->MakeTrack();
366 
367  TEveViewer *ev = gEve->GetDefaultViewer();
368  TGLViewer *gv = ev->GetGLViewer();
369  gv->SetGuideState(TGLUtil::kAxesOrigin, kTRUE, kFALSE, 0);
370 
371  gEve->Redraw3D(kTRUE);
372  gSystem->ProcessEvents();
373 
374  gv->CurrentCamera().RotateRad(-0.5, 1.4);
375  gv->RequestDraw();
376 }