ROOT   Reference Guide
Searching...
No Matches
TGeoHelix.cxx
Go to the documentation of this file.
1// @(#)root/geom:$Id$
2// Author: Andrei Gheata 28/04/04
3
4/*************************************************************************
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. * 9 * For the list of contributors see$ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12/** \class TGeoHelix
13\ingroup Shapes_classes
14
15Class representing a helix curve
16
17 A helix is a curve defined by the following equations:
18
19~~~ {.cpp}
20 x = (1/c) * COS(q*phi)
21 y = (1/c) * SIN(q*phi)
22 z = s * alfa
23~~~
24
25where:
26
27~~~ {.cpp}
28 c = 1/Rxy - curvature in XY plane
29 phi - phi angle
30 S = 2*PI*s - vertical separation between helix loops
31 q = +/- 1 - (+)=left-handed, (-)=right-handed
32~~~
33
34 In particular, a helix describes the trajectory of a charged particle in magnetic
35field. In such case, the helix is right-handed for negative particle charge.
36To define a helix, one must define:
37 - the curvature - positive defined
38 - the Z step made after one full turn of the helix
39 - the particle charge sign
40 - the initial particle position and direction (force normalization to unit)
41 - the magnetic field direction
42
43A helix provides:
44 - propagation to a given Z position (in global frame)
45 Double_t *point = TGeoHelix::PropagateToZ(Double_t z);
46 - propagation to an arbitrary plane, returning also the new point
47 - propagation in a geometry until the next crossed surface
48 - computation of the total track length along a helix
49*/
50
51#include "TMath.h"
52#include "TGeoShape.h"
53#include "TGeoMatrix.h"
54#include "TGeoHelix.h"
55
57
58////////////////////////////////////////////////////////////////////////////////
59/// Dummy constructor
60
62{
63 fC = 0.;
64 fS = 0.;
65 fStep = 0.;
66 fPhi = 0.;
67 fPointInit[0] = fPointInit[1] = fPointInit[2] = 0.;
68 fDirInit[0] = fDirInit[1] = fDirInit[2] = 0.;
69 fPoint[0] = fPoint[1] = fPoint[2] = 0.;
70 fDir[0] = fDir[1] = fDir[2] = 0.;
71 fB[0] = fB[1] = fB[2] = 0.;
72 fQ = 0;
73 fMatrix = nullptr;
77}
78
79////////////////////////////////////////////////////////////////////////////////
80/// Normal constructor
81
83{
84 SetXYcurvature(curvature);
85 SetHelixStep(hstep);
86 fQ = 0;
87 SetCharge(charge);
88 fStep = 0.;
89 fPhi = 0.;
90 fPointInit[0] = fPointInit[1] = fPointInit[2] = 0.;
91 fDirInit[0] = fDirInit[1] = fDirInit[2] = 0.;
92 fPoint[0] = fPoint[1] = fPoint[2] = 0.;
93 fDir[0] = fDir[1] = fDir[2] = 0.;
94 fB[0] = fB[1] = fB[2] = 0.;
95 fMatrix = new TGeoHMatrix();
99}
100
101////////////////////////////////////////////////////////////////////////////////
102/// Destructor
103
105{
106 if (fMatrix)
107 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{
117 return 1.E30;
119 Double_t step = TMath::Sqrt(2. * epsil / c);
120 return step;
121}
122
123////////////////////////////////////////////////////////////////////////////////
124/// Initialize coordinates of a point on the helix
125
127{
128 fPointInit[0] = x0;
129 fPointInit[1] = y0;
130 fPointInit[2] = z0;
132}
133
134////////////////////////////////////////////////////////////////////////////////
135/// Set initial point on the helix.
136
138{
139 InitPoint(point[0], point[1], point[2]);
140}
141
142////////////////////////////////////////////////////////////////////////////////
143/// Initialize particle direction (tangent on the helix in initial point)
144
145void TGeoHelix::InitDirection(Double_t dirx, Double_t diry, Double_t dirz, Bool_t is_normalized)
146{
147 fDirInit[0] = dirx;
148 fDirInit[1] = diry;
149 fDirInit[2] = dirz;
151 if (is_normalized)
152 return;
153 Double_t norm = 1. / TMath::Sqrt(dirx * dirx + diry * diry + dirz * dirz);
154 for (Int_t i = 0; i < 3; i++)
155 fDirInit[i] *= norm;
156}
157
158////////////////////////////////////////////////////////////////////////////////
159/// Initialize particle direction (tangent on the helix in initial point)
160
161void TGeoHelix::InitDirection(Double_t *dir, Bool_t is_normalized)
162{
163 InitDirection(dir[0], dir[1], dir[2], is_normalized);
164}
165
166////////////////////////////////////////////////////////////////////////////////
167/// Compute helix total curvature
168
170{
171 Double_t k = fC / (1. + fC * fC * fS * fS);
172 return k;
173}
174
175////////////////////////////////////////////////////////////////////////////////
176/// Set XY curvature: c = 1/Rxy
177
179{
180 fC = curvature;
182 if (fC < 0) {
183 Error("SetXYcurvature", "Curvature %f not valid. Must be positive.", fC);
184 return;
185 }
187 Warning("SetXYcurvature", "Curvature is zero. Helix is a straight line.");
189 }
190}
191
192////////////////////////////////////////////////////////////////////////////////
193/// Positive charge means left-handed helix.
194
196{
197 if (charge == 0) {
198 Error("ctor", "charge cannot be 0 - define it positive for a left-handed helix, negative otherwise");
199 return;
200 }
201 Int_t q = TMath::Sign(1, charge);
202 if (q == fQ)
203 return;
204 fQ = q;
206}
207
208////////////////////////////////////////////////////////////////////////////////
209/// Initialize particle direction (tangent on the helix in initial point)
210
211void TGeoHelix::SetField(Double_t bx, Double_t by, Double_t bz, Bool_t is_normalized)
212{
213 fB[0] = bx;
214 fB[1] = by;
215 fB[2] = bz;
217 if (is_normalized)
218 return;
219 Double_t norm = 1. / TMath::Sqrt(bx * bx + by * by + bz * bz);
220 for (Int_t i = 0; i < 3; i++)
221 fB[i] *= norm;
222}
223
224////////////////////////////////////////////////////////////////////////////////
225/// Set Z step of the helix on a complete turn. Positive or null.
226
228{
229 if (step < 0) {
230 Error("ctor", "Z step %f not valid. Must be positive.", step);
231 return;
232 }
234 fS = 0.5 * step / TMath::Pi();
235 if (fS < TGeoShape::Tolerance())
237}
238
239////////////////////////////////////////////////////////////////////////////////
240/// Reset current point/direction to initial values
241
243{
244 fStep = 0.;
245 memcpy(fPoint, fPointInit, 3 * sizeof(Double_t));
246 memcpy(fDir, fDirInit, 3 * sizeof(Double_t));
247}
248
249////////////////////////////////////////////////////////////////////////////////
250/// Make a step from current point along the helix and compute new point, direction and angle
251/// To reach a plane/ shape boundary, one has to:
252/// 1. Compute the safety to the plane/boundary
253/// 2. Define / update a helix according local field and particle state (position, direction, charge)
254/// 3. Compute the magnetic safety (maximum distance for which the field can be considered constant)
255/// 4. Call TGeoHelix::Step() having as argument the minimum between 1. and 3.
256/// 5. Repeat from 1. until the step to be made is small enough.
257/// 6. Add to the total step the distance along a straight line from the last point
258/// to the plane/shape boundary
259
261{
262 Int_t i;
263 fStep += step;
265 for (i = 0; i < 3; i++) {
266 fPoint[i] = fPointInit[i] + fStep * fDirInit[i];
267 fDir[i] = fDirInit[i];
268 }
269 return;
270 }
272 UpdateHelix();
273 Double_t r = 1. / fC;
274 fPhi = fStep / TMath::Sqrt(r * r + fS * fS);
275 Double_t vect[3];
276 vect[0] = r * TMath::Cos(fPhi);
277 vect[1] = -fQ * r * TMath::Sin(fPhi);
278 vect[2] = fS * fPhi;
280
281 Double_t ddb = fDirInit[0] * fB[0] + fDirInit[1] * fB[1] + fDirInit[2] * fB[2];
282 Double_t f = -TMath::Sqrt(1. - ddb * ddb);
283 vect[0] = f * TMath::Sin(fPhi);
284 vect[1] = fQ * f * TMath::Cos(fPhi);
285 vect[2] = ddb;
286 TMath::Normalize(vect);
288}
289
290////////////////////////////////////////////////////////////////////////////////
291/// Propagate initial point up to a given Z position in MARS.
292
294{
295 Double_t step = 0.;
296 Double_t snext = 1.E30;
297 Double_t dx, dy, dz;
298 Double_t ddn, pdn;
300 UpdateHelix();
301 dx = point[0] - fPoint[0];
302 dy = point[1] - fPoint[1];
303 dz = point[2] - fPoint[2];
304 pdn = dx * norm[0] + dy * norm[1] + dz * norm[2];
305 ddn = fDir[0] * norm[0] + fDir[1] * norm[1] + fDir[2] * norm[2];
307 // propagate straight line to plane
308 if ((pdn * ddn) <= 0)
309 return snext;
310 snext = pdn / ddn;
311 Step(snext);
312 return snext;
313 }
314
315 Double_t r = 1. / fC;
316 Double_t dist;
317 Double_t safety = TMath::Abs(pdn);
318 Double_t safestep = ComputeSafeStep();
319 snext = 1.E30;
320 Bool_t approaching = (ddn * pdn > 0) ? kTRUE : kFALSE;
321 if (approaching)
322 snext = pdn / ddn;
323 else if (safety > 2. * r)
324 return snext;
325 while (snext > safestep) {
326 dist = TMath::Max(safety, safestep);
327 Step(dist);
328 step += dist;
329 dx = point[0] - fPoint[0];
330 dy = point[1] - fPoint[1];
331 dz = point[2] - fPoint[2];
332 pdn = dx * norm[0] + dy * norm[1] + dz * norm[2];
333 ddn = fDir[0] * norm[0] + fDir[1] * norm[1] + fDir[2] * norm[2];
334 safety = TMath::Abs(pdn);
335 approaching = (ddn * pdn > 0) ? kTRUE : kFALSE;
336 snext = 1.E30;
337 if (approaching)
338 snext = pdn / ddn;
339 else if (safety > 2. * r) {
340 ResetStep();
341 return snext;
342 }
343 }
344 step += snext;
345 Step(snext);
346 return step;
347}
348
349////////////////////////////////////////////////////////////////////////////////
350/// Update the local helix matrix.
351
353{
355 fStep = 0.;
356 memcpy(fPoint, fPointInit, 3 * sizeof(Double_t));
357 memcpy(fDir, fDirInit, 3 * sizeof(Double_t));
358 Double_t rot[9];
359 Double_t tr[3];
360 Double_t ddb = fDirInit[0] * fB[0] + fDirInit[1] * fB[1] + fDirInit[2] * fB[2];
362 // helix is just a straight line
364 fMatrix->Clear();
365 return;
366 }
367 rot[2] = fB[0];
368 rot[5] = fB[1];
369 rot[8] = fB[2];
370 if (ddb < 0)
371 fS = -TMath::Abs(fS);
372 Double_t fy = -fQ * TMath::Sqrt(1. - ddb * ddb);
373 fy = 1. / fy;
374 rot[1] = fy * (fDirInit[0] - fB[0] * ddb);
375 rot[4] = fy * (fDirInit[1] - fB[1] * ddb);
376 rot[7] = fy * (fDirInit[2] - fB[2] * ddb);
377
378 rot[0] = rot[4] * rot[8] - rot[7] * rot[5];
379 rot[3] = rot[7] * rot[2] - rot[1] * rot[8];
380 rot[6] = rot[1] * rot[5] - rot[4] * rot[2];
381
382 tr[0] = fPointInit[0] - rot[0] / fC;
383 tr[1] = fPointInit[1] - rot[3] / fC;
384 tr[2] = fPointInit[2] - rot[6] / fC;
385
387 fMatrix->SetRotation(rot);
388}
#define f(i)
Definition RSha256.hxx:104
#define c(i)
Definition RSha256.hxx:101
double fy() const
constexpr Bool_t kFALSE
Definition RtypesCore.h:101
constexpr Bool_t kTRUE
Definition RtypesCore.h:100
#define ClassImp(name)
Definition Rtypes.h:377
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
float * q
Matrix class used for computing global transformations Should NOT be used for node definition.
Definition TGeoMatrix.h:458
void SetRotation(const Double_t *matrix)
Definition TGeoMatrix.h:516
void Clear(Option_t *option="") override
clear the data for this matrix
void SetTranslation(const Double_t *vect)
Definition TGeoMatrix.h:511
Class representing a helix curve.
Definition TGeoHelix.h:19
TGeoHMatrix * fMatrix
Definition TGeoHelix.h:31
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...
~TGeoHelix() override
Destructor.
Double_t fB[3]
Definition TGeoHelix.h:29
TGeoHelix()
Dummy constructor.
Definition TGeoHelix.cxx:61
Double_t fDir[3]
Definition TGeoHelix.h:28
Double_t fS
Definition TGeoHelix.h:22
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)
Double_t StepToPlane(Double_t *point, Double_t *norm)
Propagate initial point up to a given Z position in MARS.
Double_t GetTotalCurvature() const
Compute helix total curvature.
void SetXYcurvature(Double_t curvature)
Set XY curvature: c = 1/Rxy.
Double_t fPointInit[3]
Definition TGeoHelix.h:25
void SetHelixStep(Double_t hstep)
Set Z step of the helix on a complete turn. Positive or null.
@ kHelixNeedUpdate
Definition TGeoHelix.h:37
@ kHelixCircle
Definition TGeoHelix.h:37
@ kHelixStraight
Definition TGeoHelix.h:37
void Step(Double_t step)
Make a step from current point along the helix and compute new point, direction and angle To reach a ...
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)
Double_t fDirInit[3]
Definition TGeoHelix.h:26
Double_t fStep
Definition TGeoHelix.h:23
void SetCharge(Int_t charge)
Positive charge means left-handed helix.
Double_t fPoint[3]
Definition TGeoHelix.h:27
void InitPoint(Double_t x0, Double_t y0, Double_t z0)
Initialize coordinates of a point on the helix.
Double_t fC
Definition TGeoHelix.h:21
Double_t fPhi
Definition TGeoHelix.h:24
Int_t fQ
Definition TGeoHelix.h:30
void UpdateHelix()
Update the local helix matrix.
void ResetStep()
Reset current point/direction to initial values.
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
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
static Double_t Tolerance()
Definition TGeoShape.h:90
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition TObject.h:199
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition TObject.cxx:973
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition TObject.cxx:780
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:987
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition TMathBase.h:250
Float_t Normalize(Float_t v[3])
Normalize a vector v in place.
Definition TMath.cxx:518
T1 Sign(T1 a, T2 b)
Returns a value with the magnitude of a and the sign of b.
Definition TMathBase.h:175
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:662
Double_t Cos(Double_t)
Returns the cosine of an angle of x radians.
Definition TMath.h:594
constexpr Double_t Pi()
Definition TMath.h:37
Double_t Sin(Double_t)
Returns the sine of an angle of x radians.
Definition TMath.h:588
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:123