Line data Source code
1 : /**************************************************************************
2 : * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
3 : * *
4 : * Author: The ALICE Off-line Project. *
5 : * Contributors are mentioned in the code where appropriate. *
6 : * *
7 : * Permission to use, copy, modify and distribute this software and its *
8 : * documentation strictly for non-commercial purposes is hereby granted *
9 : * without fee, provided that the above copyright notice appears in all *
10 : * copies and that both the copyright notice and this permission notice *
11 : * appear in the supporting documentation. The authors make no claims *
12 : * about the suitability of this software for any purpose. It is *
13 : * provided "as is" without express or implied warranty. *
14 : **************************************************************************/
15 :
16 : /* $Id$ */
17 :
18 : ///////////////////////////////////////////////////////////////////
19 : // //
20 : // A straight line is coded as a point (3 Double_t) and //
21 : // 3 direction cosines //
22 : // //
23 : ///////////////////////////////////////////////////////////////////
24 :
25 : #include <Riostream.h>
26 : #include <TMath.h>
27 :
28 : #include "AliStrLine.h"
29 :
30 : using std::endl;
31 : using std::cout;
32 172 : ClassImp(AliStrLine)
33 :
34 : //________________________________________________________
35 : AliStrLine::AliStrLine() :
36 0 : TObject(),
37 0 : fWMatrix(0),
38 0 : fTpar(0)
39 0 : {
40 : // Default constructor
41 0 : for(Int_t i=0;i<3;i++) {
42 0 : fP0[i] = 0.;
43 0 : fSigma2P0[i] = 0.;
44 0 : fCd[i] = 0.;
45 : }
46 0 : SetIdPoints(65535,65535);
47 0 : }
48 :
49 : //________________________________________________________
50 : AliStrLine::AliStrLine(const Double_t *const point, const Double_t *const cd, Bool_t twopoints, UShort_t id1, UShort_t id2) :
51 672 : TObject(),
52 672 : fWMatrix(0),
53 672 : fTpar(0)
54 2016 : {
55 : // Standard constructor
56 : // if twopoints is true: point and cd are the 3D coordinates of
57 : // two points defininig the straight line
58 : // if twopoint is false: point represents the 3D coordinates of a point
59 : // belonging to the straight line and cd is the
60 : // direction in space
61 5376 : for(Int_t i=0;i<3;i++)
62 2016 : fSigma2P0[i] = 0.;
63 :
64 672 : if(twopoints)
65 464 : InitTwoPoints(point,cd);
66 : else
67 208 : InitDirection(point,cd);
68 :
69 672 : SetIdPoints(id1,id2);
70 1344 : }
71 :
72 :
73 : //________________________________________________________
74 : AliStrLine::AliStrLine(const Double_t *const point, const Double_t *const sig2point, const Double_t *const cd, Bool_t twopoints, UShort_t id1, UShort_t id2) :
75 0 : TObject(),
76 0 : fWMatrix(0),
77 0 : fTpar(0)
78 0 : {
79 : // Standard constructor
80 : // if twopoints is true: point and cd are the 3D coordinates of
81 : // two points defininig the straight line
82 : // if twopoint is false: point represents the 3D coordinates of a point
83 : // belonging to the straight line and cd is the
84 : // direction in space
85 0 : for(Int_t i=0;i<3;i++)
86 0 : fSigma2P0[i] = sig2point[i];
87 :
88 0 : if(twopoints)
89 0 : InitTwoPoints(point,cd);
90 : else
91 0 : InitDirection(point,cd);
92 :
93 0 : SetIdPoints(id1,id2);
94 0 : }
95 :
96 : //________________________________________________________
97 : AliStrLine::AliStrLine(const Double_t *const point, const Double_t *const sig2point, const Double_t *const wmat, const Double_t *const cd, Bool_t twopoints, UShort_t id1, UShort_t id2) :
98 412 : TObject(),
99 412 : fWMatrix(0),
100 412 : fTpar(0)
101 1236 : {
102 : // Standard constructor
103 : // if twopoints is true: point and cd are the 3D coordinates of
104 : // two points defininig the straight line
105 : // if twopoint is false: point represents the 3D coordinates of a point
106 : // belonging to the straight line and cd is the
107 : // direction in space
108 : Int_t k = 0;
109 824 : fWMatrix = new Double_t [6];
110 3296 : for(Int_t i=0;i<3;i++){
111 1236 : fSigma2P0[i] = sig2point[i];
112 16068 : for(Int_t j=0;j<3;j++)if(j>=i)fWMatrix[k++]=wmat[3*i+j];
113 : }
114 412 : if(twopoints)
115 232 : InitTwoPoints(point,cd);
116 : else
117 180 : InitDirection(point,cd);
118 :
119 412 : SetIdPoints(id1,id2);
120 824 : }
121 :
122 : //________________________________________________________
123 : AliStrLine::AliStrLine(const AliStrLine &source):
124 0 : TObject(source),
125 0 : fWMatrix(0),
126 0 : fTpar(source.fTpar)
127 0 : {
128 : //
129 : // copy constructor
130 : //
131 0 : for(Int_t i=0;i<3;i++){
132 0 : fP0[i]=source.fP0[i];
133 0 : fSigma2P0[i]=source.fSigma2P0[i];
134 0 : fCd[i]=source.fCd[i];
135 : }
136 0 : if(source.fWMatrix){
137 0 : fWMatrix = new Double_t [6];
138 0 : for(Int_t i=0;i<6;i++)fWMatrix[i]=source.fWMatrix[i];
139 0 : }
140 0 : for(Int_t i=0;i<2;i++) fIdPoint[i]=source.fIdPoint[i];
141 0 : }
142 :
143 : //________________________________________________________
144 : AliStrLine& AliStrLine::operator=(const AliStrLine& source)
145 : {
146 : // Assignment operator
147 0 : if(this !=&source){
148 0 : TObject::operator=(source);
149 0 : for(Int_t i=0;i<3;i++){
150 0 : fP0[i]=source.fP0[i];
151 0 : fSigma2P0[i]=source.fSigma2P0[i];
152 0 : fCd[i]=source.fCd[i];
153 : }
154 :
155 0 : delete [] fWMatrix;
156 0 : fWMatrix=0;
157 0 : if(source.fWMatrix){
158 0 : fWMatrix = new Double_t [6];
159 0 : for(Int_t i=0;i<6;i++)fWMatrix[i]=source.fWMatrix[i];
160 0 : }
161 0 : for(Int_t i=0;i<2;i++) fIdPoint[i]=source.fIdPoint[i];
162 0 : }
163 0 : return *this;
164 : }
165 :
166 : //________________________________________________________
167 : void AliStrLine::GetWMatrix(Double_t *wmat)const {
168 : // Getter for weighting matrix, as a [9] dim. array
169 744 : if(!fWMatrix)return;
170 : Int_t k = 0;
171 2976 : for(Int_t i=0;i<3;i++){
172 8928 : for(Int_t j=0;j<3;j++){
173 3348 : if(j>=i){
174 2232 : wmat[3*i+j]=fWMatrix[k++];
175 2232 : }
176 : else{
177 1116 : wmat[3*i+j]=wmat[3*j+i];
178 : }
179 : }
180 : }
181 744 : }
182 :
183 : //________________________________________________________
184 : void AliStrLine::SetWMatrix(const Double_t *wmat) {
185 : // Setter for weighting matrix, strating from a [9] dim. array
186 0 : if(fWMatrix)delete [] fWMatrix;
187 0 : fWMatrix = new Double_t [6];
188 : Int_t k = 0;
189 0 : for(Int_t i=0;i<3;i++){
190 0 : for(Int_t j=0;j<3;j++)if(j>=i)fWMatrix[k++]=wmat[3*i+j];
191 : }
192 0 : }
193 :
194 : //________________________________________________________
195 : void AliStrLine::InitDirection(const Double_t *const point, const Double_t *const cd)
196 : {
197 : // Initialization from a point and a direction
198 2168 : Double_t norm = cd[0]*cd[0]+cd[1]*cd[1]+cd[2]*cd[2];
199 :
200 1084 : if(norm) {
201 1084 : norm = TMath::Sqrt(1./norm);
202 8672 : for(Int_t i=0;i<3;++i) {
203 3252 : fP0[i]=point[i];
204 3252 : fCd[i]=cd[i]*norm;
205 : }
206 1084 : fTpar = 0.;
207 1084 : }
208 0 : else AliFatal("Null direction cosines!!!");
209 1084 : }
210 :
211 : //________________________________________________________
212 : void AliStrLine::InitTwoPoints(const Double_t *const pA, const Double_t *const pB)
213 : {
214 : // Initialization from the coordinates of two
215 : // points in the space
216 1392 : Double_t cd[3];
217 5568 : for(Int_t i=0;i<3;i++)cd[i] = pB[i]-pA[i];
218 696 : InitDirection(pA,cd);
219 696 : }
220 :
221 : //________________________________________________________
222 4872 : AliStrLine::~AliStrLine() {
223 : // destructor
224 1668 : if(fWMatrix)delete [] fWMatrix;
225 2436 : }
226 :
227 : //________________________________________________________
228 : void AliStrLine::PrintStatus() const {
229 : // Print current status
230 0 : cout <<"=======================================================\n";
231 0 : cout <<"Direction cosines: ";
232 0 : for(Int_t i=0;i<3;i++)cout <<fCd[i]<<"; ";
233 0 : cout <<endl;
234 0 : cout <<"Known point: ";
235 0 : for(Int_t i=0;i<3;i++)cout <<fP0[i]<<"; ";
236 0 : cout <<endl;
237 0 : cout <<"Error on known point: ";
238 0 : for(Int_t i=0;i<3;i++)cout <<TMath::Sqrt(fSigma2P0[i])<<"; ";
239 0 : cout <<endl;
240 0 : cout <<"Current value for the parameter: "<<fTpar<<endl;
241 0 : }
242 :
243 : //________________________________________________________
244 : Int_t AliStrLine::IsParallelTo(const AliStrLine *line) const {
245 : // returns 1 if lines are parallel, 0 if not paralel
246 : const Double_t prec=1e-14;
247 5440 : Double_t cd2[3];
248 2720 : line->GetCd(cd2);
249 :
250 2720 : Double_t vecpx=fCd[1]*cd2[2]-fCd[2]*cd2[1];
251 2720 : Double_t mod=TMath::Abs(fCd[1]*cd2[2])+TMath::Abs(fCd[2]*cd2[1]);
252 5232 : if(TMath::Abs(vecpx) > prec*mod) return 0;
253 :
254 208 : Double_t vecpy=-fCd[0]*cd2[2]+fCd[2]*cd2[0];
255 208 : mod=TMath::Abs(fCd[0]*cd2[2])+TMath::Abs(fCd[2]*cd2[0]);
256 208 : if(TMath::Abs(vecpy) > prec*mod) return 0;
257 :
258 208 : Double_t vecpz=fCd[0]*cd2[1]-fCd[1]*cd2[0];
259 208 : mod=TMath::Abs(fCd[0]*cd2[1])+TMath::Abs(fCd[1]*cd2[0]);
260 208 : if(TMath::Abs(vecpz) > prec) return 0;
261 :
262 208 : return 1;
263 2720 : }
264 : //________________________________________________________
265 : Int_t AliStrLine::Crossrphi(const AliStrLine *line)
266 : {
267 : // Cross 2 lines in the X-Y plane
268 : const Double_t prec=1e-14;
269 : const Double_t big=1e20;
270 0 : Double_t p2[3];
271 0 : Double_t cd2[3];
272 0 : line->GetP0(p2);
273 0 : line->GetCd(cd2);
274 0 : Double_t a=fCd[0];
275 0 : Double_t b=-cd2[0];
276 0 : Double_t c=p2[0]-fP0[0];
277 0 : Double_t d=fCd[1];
278 0 : Double_t e=-cd2[1];
279 0 : Double_t f=p2[1]-fP0[1];
280 0 : Double_t deno = a*e-b*d;
281 0 : Double_t mod=TMath::Abs(a*e)+TMath::Abs(b*d);
282 : Int_t retcode = 0;
283 0 : if(TMath::Abs(deno) > prec*mod) {
284 0 : fTpar = (c*e-b*f)/deno;
285 0 : }
286 : else {
287 0 : fTpar = big;
288 : retcode = -1;
289 : }
290 0 : return retcode;
291 0 : }
292 :
293 : //________________________________________________________
294 : Int_t AliStrLine::CrossPoints(AliStrLine *line, Double_t *point1, Double_t *point2){
295 : // Looks for the crossing point estimated starting from the
296 : // DCA segment
297 : const Double_t prec=1e-14;
298 1444 : Double_t p2[3];
299 722 : Double_t cd2[3];
300 722 : line->GetP0(p2);
301 722 : line->GetCd(cd2);
302 : Int_t i;
303 : Double_t k1 = 0;
304 : Double_t k2 = 0;
305 : Double_t a11 = 0;
306 5776 : for(i=0;i<3;i++){
307 2166 : k1+=(fP0[i]-p2[i])*fCd[i];
308 2166 : k2+=(fP0[i]-p2[i])*cd2[i];
309 2166 : a11+=fCd[i]*cd2[i];
310 : }
311 722 : Double_t a22 = -a11;
312 : Double_t a21 = 0;
313 : Double_t a12 = 0;
314 5776 : for(i=0;i<3;i++){
315 2166 : a21+=cd2[i]*cd2[i];
316 2166 : a12-=fCd[i]*fCd[i];
317 : }
318 722 : Double_t deno = a11*a22-a21*a12;
319 722 : Double_t mod = TMath::Abs(a11*a22)+TMath::Abs(a21*a12);
320 722 : if(TMath::Abs(deno) < prec*mod) return -1;
321 722 : fTpar = (a11*k2-a21*k1) / deno;
322 722 : Double_t par2 = (k1*a22-k2*a12) / deno;
323 722 : line->SetPar(par2);
324 722 : GetCurrentPoint(point1);
325 722 : line->GetCurrentPoint(point2);
326 : return 0;
327 722 : }
328 : //________________________________________________________________
329 : Int_t AliStrLine::Cross(AliStrLine *line, Double_t *point)
330 : {
331 :
332 : //Finds intersection between lines
333 1444 : Double_t point1[3];
334 722 : Double_t point2[3];
335 722 : Int_t retcod=CrossPoints(line,point1,point2);
336 722 : if(retcod==0){
337 5776 : for(Int_t i=0;i<3;i++)point[i]=(point1[i]+point2[i])/2.;
338 722 : return 0;
339 : }else{
340 0 : return -1;
341 : }
342 722 : }
343 :
344 : //___________________________________________________________
345 : Double_t AliStrLine::GetDCA(const AliStrLine *line) const
346 : {
347 : //Returns the distance of closest approach between two lines
348 : const Double_t prec=1e-14;
349 5440 : Double_t p2[3];
350 2720 : Double_t cd2[3];
351 2720 : line->GetP0(p2);
352 2720 : line->GetCd(cd2);
353 : Int_t i;
354 2720 : Int_t ispar=IsParallelTo(line);
355 2720 : if(ispar){
356 : Double_t dist1q=0,dist2=0,mod=0;
357 1664 : for(i=0;i<3;i++){
358 624 : dist1q+=(fP0[i]-p2[i])*(fP0[i]-p2[i]);
359 624 : dist2+=(fP0[i]-p2[i])*fCd[i];
360 624 : mod+=fCd[i]*fCd[i];
361 : }
362 208 : if(TMath::Abs(mod) > prec){
363 208 : dist2/=mod;
364 208 : return TMath::Sqrt(dist1q-dist2*dist2);
365 0 : }else{return -1;}
366 : }else{
367 2512 : Double_t perp[3];
368 2512 : perp[0]=fCd[1]*cd2[2]-fCd[2]*cd2[1];
369 2512 : perp[1]=-fCd[0]*cd2[2]+fCd[2]*cd2[0];
370 2512 : perp[2]=fCd[0]*cd2[1]-fCd[1]*cd2[0];
371 : Double_t mod=0,dist=0;
372 20096 : for(i=0;i<3;i++){
373 7536 : mod+=perp[i]*perp[i];
374 7536 : dist+=(fP0[i]-p2[i])*perp[i];
375 : }
376 2512 : if(TMath::Abs(mod) > prec) {
377 2512 : return TMath::Abs(dist/TMath::Sqrt(mod));
378 0 : } else return -1;
379 2512 : }
380 2720 : }
381 : //________________________________________________________
382 : void AliStrLine::GetCurrentPoint(Double_t *point) const {
383 : // Fills the array point with the current value on the line
384 12996 : for(Int_t i=0;i<3;i++)point[i]=fP0[i]+fCd[i]*fTpar;
385 1444 : }
386 :
387 : //________________________________________________________
388 : Double_t AliStrLine::GetDistFromPoint(const Double_t *point) const
389 : {
390 : // computes distance from point
391 416 : AliStrLine tmpline(point, fCd, kFALSE);
392 208 : return GetDCA(&tmpline);
393 208 : }
394 :
395 : //________________________________________________________
396 : Bool_t AliStrLine::GetParamAtRadius(Double_t r,Double_t &t1,Double_t &t2) const
397 : {
398 : // Input: radial distance from the origin (x=0, x=0) in the bending plane
399 : // Returns a boolean: kTRUE if the line crosses the cylinder of radius r
400 : // and axis coincident with the z axis. It returns kFALSE otherwise
401 : // Output: t1 and t2 in ascending order. The parameters of the line at
402 : // the two intersections with the cylinder
403 0 : Double_t p1= fCd[0]*fP0[0]+fCd[1]*fP0[1];
404 0 : Double_t p2=fCd[0]*fCd[0]+fCd[1]*fCd[1];
405 0 : Double_t delta=p1*p1-p2*(fP0[0]*fP0[0]+fP0[1]*fP0[1]-r*r);
406 0 : if(delta<0.){
407 0 : t1=-1000.;
408 0 : t2=t1;
409 0 : return kFALSE;
410 : }
411 0 : delta=TMath::Sqrt(delta);
412 0 : t1=(-p1-delta)/p2;
413 0 : t2=(-p1+delta)/p2;
414 :
415 0 : if(t2<t1){
416 : // use delta as a temporary buffer
417 : delta=t1;
418 0 : t1=t2;
419 0 : t2=delta;
420 0 : }
421 0 : if(TMath::AreEqualAbs(t1,t2,1.e-9))t1=t2;
422 0 : return kTRUE;
423 0 : }
|
