ROOT   Reference Guide
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
13/** \class TGeoHelix
14\ingroup Shapes_classes
15
16Class 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
26where:
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
36field. In such case, the helix is right-handed for negative particle charge.
37To 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
44A 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;
78}
79
80////////////////////////////////////////////////////////////////////////////////
81/// Normal constructor
82
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{
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
144void 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
158void 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
207void 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;
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;
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;
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
373 fMatrix->SetRotation(rot);
374
375}
#define f(i)
Definition: RSha256.hxx:104
#define c(i)
Definition: RSha256.hxx:101
const Bool_t kFALSE
Definition: RtypesCore.h:101
const Bool_t kTRUE
Definition: RtypesCore.h:100
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
ClassImp(TGeoHelix)
float * q
Definition: THbookFile.cxx:89
Matrix class used for computing global transformations Should NOT be used for node definition.
Definition: TGeoMatrix.h:421
void SetRotation(const Double_t *matrix)
Definition: TGeoMatrix.h:463
void Clear(Option_t *option="")
clear the data for this matrix
void SetTranslation(const Double_t *vect)
Definition: TGeoMatrix.h:462
Class representing a helix curve.
Definition: TGeoHelix.h:21
TGeoHMatrix * fMatrix
Definition: TGeoHelix.h:33
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
Double_t fB[3]
Definition: TGeoHelix.h:31
TGeoHelix()
Dummy constructor.
Definition: TGeoHelix.cxx:62
Double_t fDir[3]
Definition: TGeoHelix.h:30
Double_t fS
Definition: TGeoHelix.h:24
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
Double_t StepToPlane(Double_t *point, Double_t *norm)
Propagate initial point up to a given Z position in MARS.
Definition: TGeoHelix.cxx:285
Double_t GetTotalCurvature() const
Compute helix total curvature.
Definition: TGeoHelix.cxx:166
virtual ~TGeoHelix()
Destructor.
Definition: TGeoHelix.cxx:105
void SetXYcurvature(Double_t curvature)
Set XY curvature: c = 1/Rxy.
Definition: TGeoHelix.cxx:175
Double_t fPointInit[3]
Definition: TGeoHelix.h:27
void SetHelixStep(Double_t hstep)
Set Z step of the helix on a complete turn. Positive or null.
Definition: TGeoHelix.cxx:221
@ kHelixNeedUpdate
Definition: TGeoHelix.h:39
@ kHelixCircle
Definition: TGeoHelix.h:41
@ kHelixStraight
Definition: TGeoHelix.h:40
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
Double_t fDirInit[3]
Definition: TGeoHelix.h:28
Double_t fStep
Definition: TGeoHelix.h:25
void SetCharge(Int_t charge)
Positive charge means left-handed helix.
Definition: TGeoHelix.cxx:192
Double_t fPoint[3]
Definition: TGeoHelix.h:29
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 fC
Definition: TGeoHelix.h:23
Double_t fPhi
Definition: TGeoHelix.h:26
Int_t fQ
Definition: TGeoHelix.h:32
void UpdateHelix()
Update the local helix matrix.
Definition: TGeoHelix.cxx:339
void ResetStep()
Reset current point/direction to initial values.
Definition: TGeoHelix.cxx:235
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:363
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:339
static Double_t Tolerance()
Definition: TGeoShape.h:91
R__ALWAYS_INLINE Bool_t TestBit(UInt_t f) const
Definition: TObject.h:187
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:879
void SetBit(UInt_t f, Bool_t set)
Set or unset the user status bits as specified in f.
Definition: TObject.cxx:698
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:893
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
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:660
Double_t Cos(Double_t)
Returns the cosine of an angle of x radians.
Definition: TMath.h:592
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:586
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition: TMathBase.h:123
#define snext(osub1, osub2)
Definition: triangle.c:1168