Inserting TMath.h where required by the new version of ROOT
[u/mrichter/AliRoot.git] / STEER / AliStrLine.cxx
CommitLineData
edc97986 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 **************************************************************************/
090026bf 15
16/* $Id$ */
17
edc97986 18///////////////////////////////////////////////////////////////////
19// //
20// A straight line is coded as a point (3 Double_t) and //
21// 3 direction cosines //
22// //
23///////////////////////////////////////////////////////////////////
090026bf 24
edc97986 25#include <Riostream.h>
26#include <TTree.h>
090026bf 27#include <TMath.h>
28
edc97986 29#include "AliStrLine.h"
30
31ClassImp(AliStrLine)
32
33//________________________________________________________
fe12e09c 34AliStrLine::AliStrLine() :
35 TObject(),
36 fTpar(0),
37 fDebug(0)
38 {
edc97986 39 // Default constructor
40 for(Int_t i=0;i<3;i++) {
41 fP0[i] = 0.;
42 fCd[i] = 0.;
43 }
edc97986 44}
45
46//________________________________________________________
fe12e09c 47AliStrLine::AliStrLine(Double_t *point, Double_t *cd,Bool_t twopoints) :
48 TObject(),
49 fTpar(0),
50 fDebug(0)
51{
edc97986 52 // Standard constructor
24a0c65f 53 // if twopoints is true: point and cd are the 3D coordinates of
54 // two points defininig the straight line
55 // if twopoint is false: point represents the 3D coordinates of a point
56 // belonging to the straight line and cd is the
57 // direction in space
58 if(twopoints){
59 InitTwoPoints(point,cd);
60 }
61 else {
62 InitDirection(point,cd);
63 }
64}
65
2c9641ee 66//________________________________________________________
67AliStrLine::AliStrLine(Float_t *pointf, Float_t *cdf,Bool_t twopoints) :
68 TObject(),
69 fTpar(0),
70 fDebug(0)
71{
72 // Standard constructor - with float arguments
73 // if twopoints is true: point and cd are the 3D coordinates of
74 // two points defininig the straight line
75 // if twopoint is false: point represents the 3D coordinates of a point
76 // belonging to the straight line and cd is the
77 // direction in space
78 Double_t point[3];
79 Double_t cd[3];
80 for(Int_t i=0;i<3;i++){
81 point[i] = pointf[i];
82 cd[i] = cdf[i];
83 }
84 if(twopoints){
85 InitTwoPoints(point,cd);
86 }
87 else {
88 InitDirection(point,cd);
89 }
90}
24a0c65f 91
92//________________________________________________________
93void AliStrLine::InitDirection(Double_t *point, Double_t *cd){
94 // Initialization from a point and a direction
edc97986 95 Double_t norm = 0.;
96 for(Int_t i=0;i<3;i++)norm+=cd[i]*cd[i];
97 if(norm) {
98 norm = TMath::Sqrt(norm);
99 for(Int_t i=0;i<3;i++) cd[i]/=norm;
100 }
101 else {
102 Error("AliStrLine","Null direction cosines!!!");
103 }
104 SetP0(point);
105 SetCd(cd);
106 fTpar = 0.;
107 SetDebug();
108}
109
110//________________________________________________________
24a0c65f 111void AliStrLine::InitTwoPoints(Double_t *pA, Double_t *pB){
112 // Initialization from the coordinates of two
113 // points in the space
114 Double_t cd[3];
115 for(Int_t i=0;i<3;i++)cd[i] = pB[i]-pA[i];
116 InitDirection(pA,cd);
117}
118
119//________________________________________________________
edc97986 120AliStrLine::~AliStrLine() {
121 // destructor
122}
123
124//________________________________________________________
125void AliStrLine::PrintStatus() const {
126 // Print current status
127 cout <<"=======================================================\n";
128 cout <<"Direction cosines: ";
129 for(Int_t i=0;i<3;i++)cout <<fCd[i]<<"; ";
130 cout <<endl;
131 cout <<"Known point: ";
132 for(Int_t i=0;i<3;i++)cout <<fP0[i]<<"; ";
133 cout <<endl;
134 cout <<"Current value for the parameter: "<<fTpar<<endl;
135 cout <<" Debug flag: "<<fDebug<<endl;
136}
137
138//________________________________________________________
139Int_t AliStrLine::IsParallelTo(AliStrLine *line) const {
140 // returns 1 if lines are parallel, 0 if not paralel
141 Double_t cd2[3];
142 line->GetCd(cd2);
143 Double_t vecpx=fCd[1]*cd2[2]-fCd[2]*cd2[1];
144 if(vecpx!=0) return 0;
145 Double_t vecpy=-fCd[0]*cd2[2]+fCd[2]*cd2[0];
146 if(vecpy!=0) return 0;
147 Double_t vecpz=fCd[0]*cd2[1]-fCd[1]*cd2[0];
148 if(vecpz!=0) return 0;
149 return 1;
150}
151//________________________________________________________
152Int_t AliStrLine::Crossrphi(AliStrLine *line){
153 // Cross 2 lines in the X-Y plane
154 Double_t p2[3];
155 Double_t cd2[3];
156 line->GetP0(p2);
157 line->GetCd(cd2);
158 Double_t a=fCd[0];
159 Double_t b=-cd2[0];
160 Double_t c=p2[0]-fP0[0];
161 Double_t d=fCd[1];
162 Double_t e=-cd2[1];
163 Double_t f=p2[1]-fP0[1];
164 Double_t deno = a*e-b*d;
165 Int_t retcode = 0;
166 if(deno != 0.) {
167 fTpar = (c*e-b*f)/deno;
168 }
169 else {
170 retcode = -1;
171 }
172 return retcode;
173}
174
175//________________________________________________________
176Int_t AliStrLine::CrossPoints(AliStrLine *line, Double_t *point1, Double_t *point2){
177 // Looks for the crossing point estimated starting from the
178 // DCA segment
179 Double_t p2[3];
180 Double_t cd2[3];
181 line->GetP0(p2);
182 line->GetCd(cd2);
183 Int_t i;
184 Double_t k1 = 0;
185 for(i=0;i<3;i++)k1+=(fP0[i]-p2[i])*fCd[i];
186 Double_t k2 = 0;
187 for(i=0;i<3;i++)k2+=(fP0[i]-p2[i])*cd2[i];
188 Double_t a11 = 0;
189 for(i=0;i<3;i++)a11+=fCd[i]*cd2[i];
190 Double_t a22 = -a11;
191 Double_t a21 = 0;
192 for(i=0;i<3;i++)a21+=cd2[i]*cd2[i];
193 Double_t a12 = 0;
194 for(i=0;i<3;i++)a12-=fCd[i]*fCd[i];
195 Double_t deno = a11*a22-a21*a12;
196 if(deno == 0.) return -1;
197 fTpar = (a11*k2-a21*k1) / deno;
198 Double_t par2 = (k1*a22-k2*a12) / deno;
199 line->SetPar(par2);
200 GetCurrentPoint(point1);
201 line->GetCurrentPoint(point2);
202 return 0;
203}
204//________________________________________________________________
205Int_t AliStrLine::Cross(AliStrLine *line, Double_t *point){
206
207 //Finds intersection between lines
208 Double_t point1[3];
209 Double_t point2[3];
210 Int_t retcod=CrossPoints(line,point1,point2);
211 if(retcod==0){
212 for(Int_t i=0;i<3;i++)point[i]=(point1[i]+point2[i])/2.;
213 return 0;
214 }else{
215 return -1;
216 }
217}
218
219//___________________________________________________________
2c9641ee 220Double_t AliStrLine::GetDCA(AliStrLine *line) const{
edc97986 221 //Returns the distance of closest approach between two lines
222 Double_t p2[3];
223 Double_t cd2[3];
224 line->GetP0(p2);
225 line->GetCd(cd2);
226 Int_t i;
227 Int_t ispar=IsParallelTo(line);
228 if(ispar){
229 Double_t dist1q=0,dist2=0,mod=0;
230 for(i=0;i<3;i++){
231 dist1q+=(fP0[i]-p2[i])*(fP0[i]-p2[i]);
232 dist2+=(fP0[i]-p2[i])*fCd[i];
233 mod+=fCd[i]*fCd[i];
234 }
235 if(mod!=0){
236 dist2/=mod;
237 return TMath::Sqrt(dist1q-dist2*dist2);
238 }else{return -1;}
239 }else{
240 Double_t perp[3];
241 perp[0]=fCd[1]*cd2[2]-fCd[2]*cd2[1];
242 perp[1]=-fCd[0]*cd2[2]+fCd[2]*cd2[0];
243 perp[2]=fCd[0]*cd2[1]-fCd[1]*cd2[0];
244 Double_t mod=0,dist=0;
245 for(i=0;i<3;i++){
246 mod+=perp[i]*perp[i];
247 dist+=(fP0[i]-p2[i])*perp[i];
248 }
249 mod=sqrt(mod);
250 if(mod!=0){
251 dist/=mod;
252 return TMath::Abs(dist);
253 }else{return -1;}
254 }
255}
256//________________________________________________________
257void AliStrLine::GetCurrentPoint(Double_t *point) const {
258 // Fills the array point with the current value on the line
259 for(Int_t i=0;i<3;i++)point[i]=fP0[i]+fCd[i]*fTpar;
260}
2c9641ee 261
262//________________________________________________________
263Double_t AliStrLine::GetDistFromPoint(Double_t *point) const {
264 // computes distance from point
265 AliStrLine tmp(point,(Double_t *)fCd,kFALSE);
266 return this->GetDCA(&tmp);
267}