Logo ROOT   6.08/07
Reference Guide
TGeoHelix.cxx
Go to the documentation of this file.
1 // @(#)root/geom:$Id$
2 // Author: Andrei Gheata 28/04/04
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2000, 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 /** \class TGeoHelix
14 \ingroup Geometry_classes
15 
16 Class representing a helix curve
17 
18  A helix is a curve defined by the following equations:
19 
20 ~~~ {.cpp}
21  x = (1/c) * COS(q*phi)
22  y = (1/c) * SIN(q*phi)
23  z = s * alfa
24 ~~~
25 
26 where:
27 
28 ~~~ {.cpp}
29  c = 1/Rxy - curvature in XY plane
30  phi - phi angle
31  S = 2*PI*s - vertical separation between helix loops
32  q = +/- 1 - (+)=left-handed, (-)=right-handed
33 ~~~
34 
35  In particular, a helix describes the trajectory of a charged particle in magnetic
36 field. In such case, the helix is right-handed for negative particle charge.
37 To define a helix, one must define:
38  - the curvature - positive defined
39  - the Z step made after one full turn of the helix
40  - the particle charge sign
41  - the initial particle position and direction (force normalization to unit)
42  - the magnetic field direction
43 
44 A helix provides:
45  - propagation to a given Z position (in global frame)
46  Double_t *point = TGeoHelix::PropagateToZ(Double_t z);
47  - propagation to an arbitrary plane, returning also the new point
48  - propagation in a geometry until the next crossed surface
49  - computation of the total track length along a helix
50 */
51 
52 #include "TMath.h"
53 #include "TGeoShape.h"
54 #include "TGeoMatrix.h"
55 #include "TGeoHelix.h"
56 
58 
59 ////////////////////////////////////////////////////////////////////////////////
60 /// Dummy constructor
61 
63 {
64  fC = 0.;
65  fS = 0.;
66  fStep = 0.;
67  fPhi = 0.;
68  fPointInit[0] = fPointInit[1] = fPointInit[2] = 0.;
69  fDirInit[0] = fDirInit[1] = fDirInit[2] = 0.;
70  fPoint[0] = fPoint[1] = fPoint[2] = 0.;
71  fDir[0] = fDir[1] = fDir[2] = 0.;
72  fB[0] = fB[1] = fB[2] = 0.;
73  fQ = 0;
74  fMatrix = 0;
75  TObject::SetBit(kHelixNeedUpdate, kTRUE);
76  TObject::SetBit(kHelixStraight, kFALSE);
77  TObject::SetBit(kHelixCircle, kFALSE);
78 }
79 
80 ////////////////////////////////////////////////////////////////////////////////
81 /// Normal constructor
82 
83 TGeoHelix::TGeoHelix(Double_t curvature, Double_t hstep, Int_t charge)
84 {
85  SetXYcurvature(curvature);
86  SetHelixStep(hstep);
87  fQ = 0;
88  SetCharge(charge);
89  fStep = 0.;
90  fPhi = 0.;
91  fPointInit[0] = fPointInit[1] = fPointInit[2] = 0.;
92  fDirInit[0] = fDirInit[1] = fDirInit[2] = 0.;
93  fPoint[0] = fPoint[1] = fPoint[2] = 0.;
94  fDir[0] = fDir[1] = fDir[2] = 0.;
95  fB[0] = fB[1] = fB[2] = 0.;
96  fMatrix = new TGeoHMatrix();
100 }
101 
102 ////////////////////////////////////////////////////////////////////////////////
103 /// Destructor
104 
106 {
107  if (fMatrix) delete fMatrix;
108 }
109 
110 ////////////////////////////////////////////////////////////////////////////////
111 /// Compute safe linear step that can be made such that the error
112 /// between linear-helix extrapolation is less than EPSIL.
113 
115 {
116  if (TestBit(kHelixStraight) || TMath::Abs(fC)<TGeoShape::Tolerance()) return 1.E30;
118  Double_t step = TMath::Sqrt(2.*epsil/c);
119  return step;
120 }
121 
122 ////////////////////////////////////////////////////////////////////////////////
123 /// Initialize coordinates of a point on the helix
124 
126 {
127  fPointInit[0] = x0;
128  fPointInit[1] = y0;
129  fPointInit[2] = z0;
131 }
132 
133 ////////////////////////////////////////////////////////////////////////////////
134 /// Set initial point on the helix.
135 
137 {
138  InitPoint(point[0], point[1], point[2]);
139 }
140 
141 ////////////////////////////////////////////////////////////////////////////////
142 /// Initialize particle direction (tangent on the helix in initial point)
143 
144 void TGeoHelix::InitDirection(Double_t dirx, Double_t diry, Double_t dirz, Bool_t is_normalized)
145 {
146  fDirInit[0] = dirx;
147  fDirInit[1] = diry;
148  fDirInit[2] = dirz;
150  if (is_normalized) return;
151  Double_t norm = 1./TMath::Sqrt(dirx*dirx+diry*diry+dirz*dirz);
152  for (Int_t i=0; i<3; i++) fDirInit[i] *= norm;
153 }
154 
155 ////////////////////////////////////////////////////////////////////////////////
156 /// Initialize particle direction (tangent on the helix in initial point)
157 
158 void TGeoHelix::InitDirection(Double_t *dir, Bool_t is_normalized)
159 {
160  InitDirection(dir[0], dir[1], dir[2], is_normalized);
161 }
162 
163 ////////////////////////////////////////////////////////////////////////////////
164 /// Compute helix total curvature
165 
167 {
168  Double_t k = fC/(1.+fC*fC*fS*fS);
169  return k;
170 }
171 
172 ////////////////////////////////////////////////////////////////////////////////
173 /// Set XY curvature: c = 1/Rxy
174 
176 {
177  fC = curvature;
179  if (fC < 0) {
180  Error("SetXYcurvature", "Curvature %f not valid. Must be positive.", fC);
181  return;
182  }
184  Warning("SetXYcurvature", "Curvature is zero. Helix is a straight line.");
186  }
187 }
188 
189 ////////////////////////////////////////////////////////////////////////////////
190 /// Positive charge means left-handed helix.
191 
193 {
194  if (charge==0) {
195  Error("ctor", "charge cannot be 0 - define it positive for a left-handed helix, negative otherwise");
196  return;
197  }
198  Int_t q = TMath::Sign(1, charge);
199  if (q == fQ) return;
200  fQ = q;
202 }
203 
204 ////////////////////////////////////////////////////////////////////////////////
205 /// Initialize particle direction (tangent on the helix in initial point)
206 
207 void TGeoHelix::SetField(Double_t bx, Double_t by, Double_t bz, Bool_t is_normalized)
208 {
209  fB[0] = bx;
210  fB[1] = by;
211  fB[2] = bz;
213  if (is_normalized) return;
214  Double_t norm = 1./TMath::Sqrt(bx*bx+by*by+bz*bz);
215  for (Int_t i=0; i<3; i++) fB[i] *= norm;
216 }
217 
218 ////////////////////////////////////////////////////////////////////////////////
219 /// Set Z step of the helix on a complete turn. Positive or null.
220 
222 {
223  if (step < 0) {
224  Error("ctor", "Z step %f not valid. Must be positive.", step);
225  return;
226  }
228  fS = 0.5*step/TMath::Pi();
230 }
231 
232 ////////////////////////////////////////////////////////////////////////////////
233 /// Reset current point/direction to initial values
234 
236 {
237  fStep = 0.;
238  memcpy(fPoint, fPointInit, 3*sizeof(Double_t));
239  memcpy(fDir, fDirInit, 3*sizeof(Double_t));
240 }
241 
242 ////////////////////////////////////////////////////////////////////////////////
243 /// Make a step from current point along the helix and compute new point, direction and angle
244 /// To reach a plane/ shape boundary, one has to:
245 /// 1. Compute the safety to the plane/boundary
246 /// 2. Define / update a helix according local field and particle state (position, direction, charge)
247 /// 3. Compute the magnetic safety (maximum distance for which the field can be considered constant)
248 /// 4. Call TGeoHelix::Step() having as argument the minimum between 1. and 3.
249 /// 5. Repeat from 1. until the step to be made is small enough.
250 /// 6. Add to the total step the distance along a straight line from the last point
251 /// to the plane/shape boundary
252 
254 {
255  Int_t i;
256  fStep += step;
258  for (i=0; i<3; i++) {
259  fPoint[i] = fPointInit[i]+fStep*fDirInit[i];
260  fDir[i] = fDirInit[i];
261  }
262  return;
263  }
265  Double_t r = 1./fC;
266  fPhi = fStep/TMath::Sqrt(r*r+fS*fS);
267  Double_t vect[3];
268  vect[0] = r * TMath::Cos(fPhi);
269  vect[1] = -fQ * r * TMath::Sin(fPhi);
270  vect[2] = fS * fPhi;
271  fMatrix->LocalToMaster(vect, fPoint);
272 
273  Double_t ddb = fDirInit[0]*fB[0]+fDirInit[1]*fB[1]+fDirInit[2]*fB[2];
274  Double_t f = -TMath::Sqrt(1.-ddb*ddb);
275  vect[0] = f*TMath::Sin(fPhi);
276  vect[1] = fQ*f*TMath::Cos(fPhi);
277  vect[2] = ddb;
278  TMath::Normalize(vect);
280 }
281 
282 ////////////////////////////////////////////////////////////////////////////////
283 /// Propagate initial point up to a given Z position in MARS.
284 
286 {
287  Double_t step = 0.;
288  Double_t snext = 1.E30;
289  Double_t dx, dy, dz;
290  Double_t ddn, pdn;
292  dx = point[0] - fPoint[0];
293  dy = point[1] - fPoint[1];
294  dz = point[2] - fPoint[2];
295  pdn = dx*norm[0]+dy*norm[1]+dz*norm[2];
296  ddn = fDir[0]*norm[0]+fDir[1]*norm[1]+fDir[2]*norm[2];
298  // propagate straight line to plane
299  if ((pdn*ddn) <= 0) return snext;
300  snext = pdn/ddn;
301  Step(snext);
302  return snext;
303  }
304 
305  Double_t r = 1./fC;
306  Double_t dist;
307  Double_t safety = TMath::Abs(pdn);
308  Double_t safestep = ComputeSafeStep();
309  snext = 1.E30;
310  Bool_t approaching = (ddn*pdn>0)?kTRUE:kFALSE;
311  if (approaching) snext = pdn/ddn;
312  else if (safety > 2.*r) return snext;
313  while (snext > safestep) {
314  dist = TMath::Max(safety, safestep);
315  Step(dist);
316  step += dist;
317  dx = point[0] - fPoint[0];
318  dy = point[1] - fPoint[1];
319  dz = point[2] - fPoint[2];
320  pdn = dx*norm[0]+dy*norm[1]+dz*norm[2];
321  ddn = fDir[0]*norm[0]+fDir[1]*norm[1]+fDir[2]*norm[2];
322  safety = TMath::Abs(pdn);
323  approaching = (ddn*pdn>0)?kTRUE:kFALSE;
324  snext = 1.E30;
325  if (approaching) snext = pdn/ddn;
326  else if (safety > 2.*r) {
327  ResetStep();
328  return snext;
329  }
330  }
331  step += snext;
332  Step(snext);
333  return step;
334 }
335 
336 ////////////////////////////////////////////////////////////////////////////////
337 /// Update the local helix matrix.
338 
340 {
342  fStep = 0.;
343  memcpy(fPoint, fPointInit, 3*sizeof(Double_t));
344  memcpy(fDir, fDirInit, 3*sizeof(Double_t));
345  Double_t rot[9];
346  Double_t tr[3];
347  Double_t ddb = fDirInit[0]*fB[0]+fDirInit[1]*fB[1]+fDirInit[2]*fB[2];
349  // helix is just a straight line
351  fMatrix->Clear();
352  return;
353  }
354  rot[2] = fB[0];
355  rot[5] = fB[1];
356  rot[8] = fB[2];
357  if (ddb < 0) fS = -TMath::Abs(fS);
358  Double_t fy = - fQ*TMath::Sqrt(1.-ddb*ddb);
359  fy = 1./fy;
360  rot[1] = fy*(fDirInit[0]-fB[0]*ddb);
361  rot[4] = fy*(fDirInit[1]-fB[1]*ddb);
362  rot[7] = fy*(fDirInit[2]-fB[2]*ddb);
363 
364  rot[0] = rot[4]*rot[8] - rot[7]*rot[5];
365  rot[3] = rot[7]*rot[2] - rot[1]*rot[8];
366  rot[6] = rot[1]*rot[5] - rot[4]*rot[2];
367 
368  tr[0] = fPointInit[0] - rot[0]/fC;
369  tr[1] = fPointInit[1] - rot[3]/fC;
370  tr[2] = fPointInit[2] - rot[6]/fC;
371 
372  fMatrix->SetTranslation(tr);
373  fMatrix->SetRotation(rot);
374 
375 }
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
#define snext(osub1, osub2)
Definition: triangle.c:1167
Double_t GetTotalCurvature() const
Compute helix total curvature.
Definition: TGeoHelix.cxx:166
Float_t Normalize(Float_t v[3])
Normalize a vector v in place.
Definition: TMath.cxx:498
T1 Sign(T1 a, T2 b)
Definition: TMathBase.h:155
return c
void SetCharge(Int_t charge)
Positive charge means left-handed helix.
Definition: TGeoHelix.cxx:192
Bool_t TestBit(UInt_t f) const
Definition: TObject.h:157
void SetTranslation(const Double_t *vect)
Definition: TGeoMatrix.h:448
void UpdateHelix()
Update the local helix matrix.
Definition: TGeoHelix.cxx:339
Matrix class used for computing global transformations Should NOT be used for node definition...
Definition: TGeoMatrix.h:410
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
void SetHelixStep(Double_t hstep)
Set Z step of the helix on a complete turn. Positive or null.
Definition: TGeoHelix.cxx:221
const Bool_t kFALSE
Definition: Rtypes.h:92
TGeoHelix()
Dummy constructor.
Definition: TGeoHelix.cxx:62
void Step(Double_t step)
Make a step from current point along the helix and compute new point, direction and angle To reach a ...
Definition: TGeoHelix.cxx:253
void SetField(Double_t bx, Double_t by, Double_t bz, Bool_t is_normalized=kTRUE)
Initialize particle direction (tangent on the helix in initial point)
Definition: TGeoHelix.cxx:207
Short_t Abs(Short_t d)
Definition: TMathBase.h:110
Double_t StepToPlane(Double_t *point, Double_t *norm)
Propagate initial point up to a given Z position in MARS.
Definition: TGeoHelix.cxx:285
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:739
Double_t fB[3]
Definition: TGeoHelix.h:33
static Double_t Tolerance()
Definition: TGeoShape.h:93
Double_t fPhi
Definition: TGeoHelix.h:28
Double_t fPoint[3]
Definition: TGeoHelix.h:31
virtual ~TGeoHelix()
Destructor.
Definition: TGeoHelix.cxx:105
Int_t fQ
Definition: TGeoHelix.h:34
Double_t fC
Definition: TGeoHelix.h:25
virtual void LocalToMasterVect(const Double_t *local, Double_t *master) const
convert a vector by multiplying its column vector (x, y, z, 1) to matrix inverse
Definition: TGeoMatrix.cxx:413
Double_t fStep
Definition: TGeoHelix.h:27
void ResetStep()
Reset current point/direction to initial values.
Definition: TGeoHelix.cxx:235
TRandom2 r(17)
void Clear(Option_t *option="")
clear the data for this matrix
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:925
Double_t fS
Definition: TGeoHelix.h:26
virtual void LocalToMaster(const Double_t *local, Double_t *master) const
convert a point by multiplying its column vector (x, y, z, 1) to matrix inverse
Definition: TGeoMatrix.cxx:389
Double_t Cos(Double_t)
Definition: TMath.h:424
void SetXYcurvature(Double_t curvature)
Set XY curvature: c = 1/Rxy.
Definition: TGeoHelix.cxx:175
Double_t Pi()
Definition: TMath.h:44
Double_t fDir[3]
Definition: TGeoHelix.h:32
void InitDirection(Double_t dirx, Double_t diry, Double_t dirz, Bool_t is_normalized=kTRUE)
Initialize particle direction (tangent on the helix in initial point)
Definition: TGeoHelix.cxx:144
#define ClassImp(name)
Definition: Rtypes.h:279
double f(double x)
Class representing a helix curve.
Definition: TGeoHelix.h:22
double Double_t
Definition: RtypesCore.h:55
Double_t fPointInit[3]
Definition: TGeoHelix.h:29
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:202
void SetRotation(const Double_t *matrix)
Definition: TGeoMatrix.h:449
TGeoHMatrix * fMatrix
Definition: TGeoHelix.h:35
Double_t Sin(Double_t)
Definition: TMath.h:421
Double_t ComputeSafeStep(Double_t epsil=1E-6) const
Compute safe linear step that can be made such that the error between linear-helix extrapolation is l...
Definition: TGeoHelix.cxx:114
void InitPoint(Double_t x0, Double_t y0, Double_t z0)
Initialize coordinates of a point on the helix.
Definition: TGeoHelix.cxx:125
Double_t Sqrt(Double_t x)
Definition: TMath.h:464
const Bool_t kTRUE
Definition: Rtypes.h:91
float * q
Definition: THbookFile.cxx:87
double norm(double *x, double *p)
Definition: unuranDistr.cxx:40
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:911
Double_t fDirInit[3]
Definition: TGeoHelix.h:30