Logo ROOT  
Reference Guide
 
Loading...
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/*************************************************************************
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/** \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
56
57////////////////////////////////////////////////////////////////////////////////
58/// Dummy constructor
59
61{
62 fC = 0.;
63 fS = 0.;
64 fStep = 0.;
65 fPhi = 0.;
66 fPointInit[0] = fPointInit[1] = fPointInit[2] = 0.;
67 fDirInit[0] = fDirInit[1] = fDirInit[2] = 0.;
68 fPoint[0] = fPoint[1] = fPoint[2] = 0.;
69 fDir[0] = fDir[1] = fDir[2] = 0.;
70 fB[0] = fB[1] = fB[2] = 0.;
71 fQ = 0;
72 fMatrix = nullptr;
76}
77
78////////////////////////////////////////////////////////////////////////////////
79/// Normal constructor
80
82{
85 fQ = 0;
87 fStep = 0.;
88 fPhi = 0.;
89 fPointInit[0] = fPointInit[1] = fPointInit[2] = 0.;
90 fDirInit[0] = fDirInit[1] = fDirInit[2] = 0.;
91 fPoint[0] = fPoint[1] = fPoint[2] = 0.;
92 fDir[0] = fDir[1] = fDir[2] = 0.;
93 fB[0] = fB[1] = fB[2] = 0.;
94 fMatrix = new TGeoHMatrix();
98}
99
100////////////////////////////////////////////////////////////////////////////////
101/// Destructor
102
104{
105 if (fMatrix)
106 delete fMatrix;
107}
108
109////////////////////////////////////////////////////////////////////////////////
110/// Compute safe linear step that can be made such that the error
111/// between linear-helix extrapolation is less than EPSIL.
112
114{
116 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
145{
146 fDirInit[0] = dirx;
147 fDirInit[1] = diry;
148 fDirInit[2] = dirz;
150 if (is_normalized)
151 return;
152 Double_t norm = 1. / TMath::Sqrt(dirx * dirx + diry * diry + dirz * dirz);
153 for (Int_t i = 0; i < 3; i++)
154 fDirInit[i] *= norm;
155}
156
157////////////////////////////////////////////////////////////////////////////////
158/// Initialize particle direction (tangent on the helix in initial point)
159
161{
162 InitDirection(dir[0], dir[1], dir[2], is_normalized);
163}
164
165////////////////////////////////////////////////////////////////////////////////
166/// Compute helix total curvature
167
169{
170 Double_t k = fC / (1. + fC * fC * fS * fS);
171 return k;
172}
173
174////////////////////////////////////////////////////////////////////////////////
175/// Set XY curvature: c = 1/Rxy
176
178{
179 fC = curvature;
181 if (fC < 0) {
182 Error("SetXYcurvature", "Curvature %f not valid. Must be positive.", fC);
183 return;
184 }
186 Warning("SetXYcurvature", "Curvature is zero. Helix is a straight line.");
188 }
189}
190
191////////////////////////////////////////////////////////////////////////////////
192/// Positive charge means left-handed helix.
193
195{
196 if (charge == 0) {
197 Error("ctor", "charge cannot be 0 - define it positive for a left-handed helix, negative otherwise");
198 return;
199 }
201 if (q == fQ)
202 return;
203 fQ = q;
205}
206
207////////////////////////////////////////////////////////////////////////////////
208/// Initialize particle direction (tangent on the helix in initial point)
209
211{
212 fB[0] = bx;
213 fB[1] = by;
214 fB[2] = bz;
216 if (is_normalized)
217 return;
218 Double_t norm = 1. / TMath::Sqrt(bx * bx + by * by + bz * bz);
219 for (Int_t i = 0; i < 3; i++)
220 fB[i] *= norm;
221}
222
223////////////////////////////////////////////////////////////////////////////////
224/// Set Z step of the helix on a complete turn. Positive or null.
225
227{
228 if (step < 0) {
229 Error("ctor", "Z step %f not valid. Must be positive.", step);
230 return;
231 }
233 fS = 0.5 * step / TMath::Pi();
234 if (fS < TGeoShape::Tolerance())
236}
237
238////////////////////////////////////////////////////////////////////////////////
239/// Reset current point/direction to initial values
240
242{
243 fStep = 0.;
244 memcpy(fPoint, fPointInit, 3 * sizeof(Double_t));
245 memcpy(fDir, fDirInit, 3 * sizeof(Double_t));
246}
247
248////////////////////////////////////////////////////////////////////////////////
249/// Make a step from current point along the helix and compute new point, direction and angle
250/// To reach a plane/ shape boundary, one has to:
251/// 1. Compute the safety to the plane/boundary
252/// 2. Define / update a helix according local field and particle state (position, direction, charge)
253/// 3. Compute the magnetic safety (maximum distance for which the field can be considered constant)
254/// 4. Call TGeoHelix::Step() having as argument the minimum between 1. and 3.
255/// 5. Repeat from 1. until the step to be made is small enough.
256/// 6. Add to the total step the distance along a straight line from the last point
257/// to the plane/shape boundary
258
260{
261 Int_t i;
262 fStep += step;
264 for (i = 0; i < 3; i++) {
265 fPoint[i] = fPointInit[i] + fStep * fDirInit[i];
266 fDir[i] = fDirInit[i];
267 }
268 return;
269 }
271 UpdateHelix();
272 Double_t r = 1. / fC;
273 fPhi = fStep / TMath::Sqrt(r * r + fS * fS);
274 Double_t vect[3];
275 vect[0] = r * TMath::Cos(fPhi);
276 vect[1] = -fQ * r * TMath::Sin(fPhi);
277 vect[2] = fS * fPhi;
279
280 Double_t ddb = fDirInit[0] * fB[0] + fDirInit[1] * fB[1] + fDirInit[2] * fB[2];
281 Double_t f = -TMath::Sqrt(1. - ddb * ddb);
282 vect[0] = f * TMath::Sin(fPhi);
283 vect[1] = fQ * f * TMath::Cos(fPhi);
284 vect[2] = ddb;
287}
288
289////////////////////////////////////////////////////////////////////////////////
290/// Propagate initial point up to a given Z position in MARS.
291
293{
294 Double_t step = 0.;
295 Double_t snext = 1.E30;
296 Double_t dx, dy, dz;
299 UpdateHelix();
300 dx = point[0] - fPoint[0];
301 dy = point[1] - fPoint[1];
302 dz = point[2] - fPoint[2];
303 pdn = dx * norm[0] + dy * norm[1] + dz * norm[2];
304 ddn = fDir[0] * norm[0] + fDir[1] * norm[1] + fDir[2] * norm[2];
306 // propagate straight line to plane
307 if ((pdn * ddn) <= 0)
308 return snext;
309 snext = pdn / ddn;
310 Step(snext);
311 return snext;
312 }
313
314 Double_t r = 1. / fC;
315 Double_t dist;
318 snext = 1.E30;
319 Bool_t approaching = (ddn * pdn > 0) ? kTRUE : kFALSE;
320 if (approaching)
321 snext = pdn / ddn;
322 else if (safety > 2. * r)
323 return snext;
324 while (snext > safestep) {
325 dist = TMath::Max(safety, safestep);
326 Step(dist);
327 step += dist;
328 dx = point[0] - fPoint[0];
329 dy = point[1] - fPoint[1];
330 dz = point[2] - fPoint[2];
331 pdn = dx * norm[0] + dy * norm[1] + dz * norm[2];
332 ddn = fDir[0] * norm[0] + fDir[1] * norm[1] + fDir[2] * norm[2];
334 approaching = (ddn * pdn > 0) ? kTRUE : kFALSE;
335 snext = 1.E30;
336 if (approaching)
337 snext = pdn / ddn;
338 else if (safety > 2. * r) {
339 ResetStep();
340 return snext;
341 }
342 }
343 step += snext;
344 Step(snext);
345 return step;
346}
347
348////////////////////////////////////////////////////////////////////////////////
349/// Update the local helix matrix.
350
352{
354 fStep = 0.;
355 memcpy(fPoint, fPointInit, 3 * sizeof(Double_t));
356 memcpy(fDir, fDirInit, 3 * sizeof(Double_t));
357 Double_t rot[9];
358 Double_t tr[3];
359 Double_t ddb = fDirInit[0] * fB[0] + fDirInit[1] * fB[1] + fDirInit[2] * fB[2];
361 // helix is just a straight line
363 fMatrix->Clear();
364 return;
365 }
366 rot[2] = fB[0];
367 rot[5] = fB[1];
368 rot[8] = fB[2];
369 if (ddb < 0)
370 fS = -TMath::Abs(fS);
371 Double_t fy = -fQ * TMath::Sqrt(1. - ddb * ddb);
372 fy = 1. / fy;
373 rot[1] = fy * (fDirInit[0] - fB[0] * ddb);
374 rot[4] = fy * (fDirInit[1] - fB[1] * ddb);
375 rot[7] = fy * (fDirInit[2] - fB[2] * ddb);
376
377 rot[0] = rot[4] * rot[8] - rot[7] * rot[5];
378 rot[3] = rot[7] * rot[2] - rot[1] * rot[8];
379 rot[6] = rot[1] * rot[5] - rot[4] * rot[2];
380
381 tr[0] = fPointInit[0] - rot[0] / fC;
382 tr[1] = fPointInit[1] - rot[3] / fC;
383 tr[2] = fPointInit[2] - rot[6] / fC;
384
387}
#define f(i)
Definition RSha256.hxx:104
#define c(i)
Definition RSha256.hxx:101
constexpr Bool_t kFALSE
Definition RtypesCore.h:108
constexpr Bool_t kTRUE
Definition RtypesCore.h:107
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
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
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:60
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:97
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition TObject.h:202
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition TObject.cxx:1057
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition TObject.cxx:864
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:1071
Short_t Max(Short_t a, Short_t b)
Returns the largest of a and b.
Definition TMathBase.h:251
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:176
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:673
Double_t Cos(Double_t)
Returns the cosine of an angle of x radians.
Definition TMath.h:605
constexpr Double_t Pi()
Definition TMath.h:40
Double_t Sin(Double_t)
Returns the sine of an angle of x radians.
Definition TMath.h:599
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:124