Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
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 
57 ClassImp(TGeoHelix);
58 
59 ////////////////////////////////////////////////////////////////////////////////
60 /// Dummy constructor
61 
62 TGeoHelix::TGeoHelix()
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();
97  TObject::SetBit(kHelixNeedUpdate, kTRUE);
98  TObject::SetBit(kHelixStraight, kFALSE);
99  TObject::SetBit(kHelixCircle, kFALSE);
100 }
101 
102 ////////////////////////////////////////////////////////////////////////////////
103 /// Destructor
104 
105 TGeoHelix::~TGeoHelix()
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 
114 Double_t TGeoHelix::ComputeSafeStep(Double_t epsil) const
115 {
116  if (TestBit(kHelixStraight) || TMath::Abs(fC)<TGeoShape::Tolerance()) return 1.E30;
117  Double_t c = GetTotalCurvature();
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 
125 void TGeoHelix::InitPoint(Double_t x0, Double_t y0, Double_t z0)
126 {
127  fPointInit[0] = x0;
128  fPointInit[1] = y0;
129  fPointInit[2] = z0;
130  TObject::SetBit(kHelixNeedUpdate, kTRUE);
131 }
132 
133 ////////////////////////////////////////////////////////////////////////////////
134 /// Set initial point on the helix.
135 
136 void TGeoHelix::InitPoint (Double_t *point)
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;
149  TObject::SetBit(kHelixNeedUpdate, kTRUE);
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 
166 Double_t TGeoHelix::GetTotalCurvature() const
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 
175 void TGeoHelix::SetXYcurvature(Double_t curvature)
176 {
177  fC = curvature;
178  TObject::SetBit(kHelixNeedUpdate, kTRUE);
179  if (fC < 0) {
180  Error("SetXYcurvature", "Curvature %f not valid. Must be positive.", fC);
181  return;
182  }
183  if (TMath::Abs(fC) < TGeoShape::Tolerance()) {
184  Warning("SetXYcurvature", "Curvature is zero. Helix is a straight line.");
185  TObject::SetBit(kHelixStraight, kTRUE);
186  }
187 }
188 
189 ////////////////////////////////////////////////////////////////////////////////
190 /// Positive charge means left-handed helix.
191 
192 void TGeoHelix::SetCharge(Int_t charge)
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;
201  TObject::SetBit(kHelixNeedUpdate, kTRUE);
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;
212  TObject::SetBit(kHelixNeedUpdate, kTRUE);
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 
221 void TGeoHelix::SetHelixStep(Double_t step)
222 {
223  if (step < 0) {
224  Error("ctor", "Z step %f not valid. Must be positive.", step);
225  return;
226  }
227  TObject::SetBit(kHelixNeedUpdate, kTRUE);
228  fS = 0.5*step/TMath::Pi();
229  if (fS < TGeoShape::Tolerance()) TObject::SetBit(kHelixCircle, kTRUE);
230 }
231 
232 ////////////////////////////////////////////////////////////////////////////////
233 /// Reset current point/direction to initial values
234 
235 void TGeoHelix::ResetStep()
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 
253 void TGeoHelix::Step(Double_t step)
254 {
255  Int_t i;
256  fStep += step;
257  if (TObject::TestBit(kHelixStraight)) {
258  for (i=0; i<3; i++) {
259  fPoint[i] = fPointInit[i]+fStep*fDirInit[i];
260  fDir[i] = fDirInit[i];
261  }
262  return;
263  }
264  if (TObject::TestBit(kHelixNeedUpdate)) UpdateHelix();
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);
279  fMatrix->LocalToMasterVect(vect, fDir);
280 }
281 
282 ////////////////////////////////////////////////////////////////////////////////
283 /// Propagate initial point up to a given Z position in MARS.
284 
285 Double_t TGeoHelix::StepToPlane(Double_t *point, Double_t *norm)
286 {
287  Double_t step = 0.;
288  Double_t snext = 1.E30;
289  Double_t dx, dy, dz;
290  Double_t ddn, pdn;
291  if (TObject::TestBit(kHelixNeedUpdate)) UpdateHelix();
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];
297  if (TObject::TestBit(kHelixStraight)) {
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 
339 void TGeoHelix::UpdateHelix()
340 {
341  TObject::SetBit(kHelixNeedUpdate, kFALSE);
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];
348  if ((1.-TMath::Abs(ddb))<TGeoShape::Tolerance() || TMath::Abs(fC)<TGeoShape::Tolerance()) {
349  // helix is just a straight line
350  TObject::SetBit(kHelixStraight, kTRUE);
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 }