- Three classes by MinJung Kweon AliHFEpriVtx, AliHFEsecVtx and AliHFEmcQA for primar...
[u/mrichter/AliRoot.git] / STEER / AliStrLine.cxx
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 <TTree.h>
27 #include <TMath.h>
28
29 #include "AliLog.h"
30 #include "AliStrLine.h"
31
32 ClassImp(AliStrLine)
33
34 //________________________________________________________
35 AliStrLine::AliStrLine() :
36   TObject(),
37   fWMatrix(0),
38   fTpar(0)
39  {
40   // Default constructor
41   for(Int_t i=0;i<3;i++) {
42     fP0[i] = 0.;
43     fSigma2P0[i] = 0.;
44     fCd[i] = 0.;
45   }
46   SetIdPoints(65535,65535);
47 }
48
49 //________________________________________________________
50 AliStrLine::AliStrLine(Double_t *point, Double_t *cd,Bool_t twopoints, UShort_t id1, UShort_t id2) :
51   TObject(),
52   fWMatrix(0),
53   fTpar(0)
54 {
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   for(Int_t i=0;i<3;i++){ 
62     fSigma2P0[i] = 0.;
63   }
64   if(twopoints){
65     InitTwoPoints(point,cd);
66   }
67   else {
68     InitDirection(point,cd);
69   }
70   SetIdPoints(id1,id2);
71 }
72
73 //________________________________________________________
74 AliStrLine::AliStrLine(Float_t *pointf, Float_t *cdf,Bool_t twopoints, UShort_t id1, UShort_t id2) :
75   TObject(),
76   fWMatrix(0),
77   fTpar(0)
78 {
79   // Standard constructor - with float arguments
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   Double_t point[3];
86   Double_t cd[3];
87   for(Int_t i=0;i<3;i++){
88     point[i] = pointf[i];
89     cd[i] = cdf[i];
90     fSigma2P0[i] = 0.;
91   }
92   if(twopoints){
93     InitTwoPoints(point,cd);
94   }
95   else {
96     InitDirection(point,cd);
97   }
98   SetIdPoints(id1,id2);
99 }
100
101 //________________________________________________________
102 AliStrLine::AliStrLine(Double_t *point, Double_t *sig2point, Double_t *cd,Bool_t twopoints, UShort_t id1, UShort_t id2) :
103   TObject(),
104   fWMatrix(0),
105   fTpar(0)
106 {
107   // Standard constructor
108   // if twopoints is true:  point and cd are the 3D coordinates of
109   //                        two points defininig the straight line
110   // if twopoint is false: point represents the 3D coordinates of a point
111   //                       belonging to the straight line and cd is the
112   //                       direction in space
113   for(Int_t i=0;i<3;i++){ 
114     fSigma2P0[i] = sig2point[i];
115   }
116   if(twopoints){
117     InitTwoPoints(point,cd);
118   }
119   else {
120     InitDirection(point,cd);
121   }
122   SetIdPoints(id1,id2);
123 }
124
125 //________________________________________________________
126 AliStrLine::AliStrLine(Float_t *pointf, Float_t *sig2point, Float_t *cdf,Bool_t twopoints, UShort_t id1, UShort_t id2) :
127   TObject(),
128   fWMatrix(0),
129   fTpar(0)
130 {
131   // Standard constructor - with float arguments
132   // if twopoints is true:  point and cd are the 3D coordinates of
133   //                        two points defininig the straight line
134   // if twopoint is false: point represents the 3D coordinates of a point
135   //                       belonging to the straight line and cd is the
136   //                       direction in space
137   Double_t point[3];
138   Double_t cd[3];
139   for(Int_t i=0;i<3;i++){
140     point[i] = pointf[i];
141     cd[i] = cdf[i];
142     fSigma2P0[i] = sig2point[i];
143   }
144   if(twopoints){
145     InitTwoPoints(point,cd);
146   }
147   else {
148     InitDirection(point,cd);
149   }
150   SetIdPoints(id1,id2);
151 }
152 //________________________________________________________
153 AliStrLine::AliStrLine(Double_t *point, Double_t *sig2point, Double_t *wmat, Double_t *cd,Bool_t twopoints, UShort_t id1, UShort_t id2) :
154   TObject(),
155   fWMatrix(0),
156   fTpar(0)
157 {
158   // Standard constructor
159   // if twopoints is true:  point and cd are the 3D coordinates of
160   //                        two points defininig the straight line
161   // if twopoint is false: point represents the 3D coordinates of a point
162   //                       belonging to the straight line and cd is the
163   //                       direction in space
164   Int_t k = 0;
165   fWMatrix = new Double_t [6];
166   for(Int_t i=0;i<3;i++){ 
167     fSigma2P0[i] = sig2point[i];
168     for(Int_t j=0;j<3;j++)if(j>=i)fWMatrix[k++]=wmat[3*i+j];
169   }
170   if(twopoints){
171     InitTwoPoints(point,cd);
172   }
173   else {
174     InitDirection(point,cd);
175   }
176   SetIdPoints(id1,id2);
177 }
178
179 //________________________________________________________
180 AliStrLine::AliStrLine(Float_t *pointf, Float_t *sig2point, Float_t *wmat, Float_t *cdf,Bool_t twopoints, UShort_t id1, UShort_t id2) :
181   TObject(),
182   fWMatrix(0),
183   fTpar(0)
184 {
185   // Standard constructor - with float arguments
186   // if twopoints is true:  point and cd are the 3D coordinates of
187   //                        two points defininig the straight line
188   // if twopoint is false: point represents the 3D coordinates of a point
189   //                       belonging to the straight line and cd is the
190   //                       direction in space
191   Double_t point[3];
192   Double_t cd[3];
193   fWMatrix = new Double_t [6];
194   Int_t k = 0;
195   for(Int_t i=0;i<3;i++){
196     point[i] = pointf[i];
197     cd[i] = cdf[i];
198     fSigma2P0[i] = sig2point[i];
199     for(Int_t j=0;j<3;j++)if(j>=i)fWMatrix[k++]=wmat[3*i+j];
200   }
201   if(twopoints){
202     InitTwoPoints(point,cd);
203   }
204   else {
205     InitDirection(point,cd);
206   }
207   SetIdPoints(id1,id2);
208 }
209
210 //________________________________________________________
211 AliStrLine::AliStrLine(const AliStrLine &source):TObject(source),
212 fWMatrix(0),
213 fTpar(source.fTpar){
214   // copy constructor
215   for(Int_t i=0;i<3;i++){
216     fP0[i]=source.fP0[i];
217     fSigma2P0[i]=source.fSigma2P0[i];
218     fCd[i]=source.fCd[i];
219   }
220   if(source.fWMatrix){
221     fWMatrix = new Double_t [6];
222     for(Int_t i=0;i<6;i++)fWMatrix[i]=source.fWMatrix[i];
223   }
224   for(Int_t i=0;i<2;i++) fIdPoint[i]=source.fIdPoint[i];
225 }
226
227 //________________________________________________________
228 AliStrLine& AliStrLine::operator=(const AliStrLine& source){
229   // Assignment operator
230   if(this !=&source){
231     this->~AliStrLine();
232     new(this)AliStrLine(source);
233   }
234   return *this;
235 }
236
237 //________________________________________________________
238 void AliStrLine::GetWMatrix(Double_t *wmat)const {
239 // Getter for weighting matrix, as a [9] dim. array
240   if(!fWMatrix)return;
241   Int_t k = 0;
242   for(Int_t i=0;i<3;i++){
243     for(Int_t j=0;j<3;j++){
244       if(j>=i){
245         wmat[3*i+j]=fWMatrix[k++];
246       }
247       else{
248         wmat[3*i+j]=wmat[3*j+i];
249       }
250     }
251   }
252
253
254 //________________________________________________________
255 void AliStrLine::SetWMatrix(const Double_t *wmat) {
256 // Setter for weighting matrix, strating from a [9] dim. array
257   if(fWMatrix)delete [] fWMatrix;
258   fWMatrix = new Double_t [6];
259   Int_t k = 0;
260   for(Int_t i=0;i<3;i++){
261     for(Int_t j=0;j<3;j++)if(j>=i)fWMatrix[k++]=wmat[3*i+j];
262   }
263 }
264
265 //________________________________________________________
266 void AliStrLine::InitDirection(Double_t *point, Double_t *cd){
267   // Initialization from a point and a direction
268   Double_t norm = 0.;
269   for(Int_t i=0;i<3;i++)norm+=cd[i]*cd[i];
270   if(norm) {
271     norm = TMath::Sqrt(norm);
272     for(Int_t i=0;i<3;i++) cd[i]/=norm;
273   }
274   else {
275     Error("AliStrLine","Null direction cosines!!!");
276   }
277   SetP0(point);
278   SetCd(cd);
279   fTpar = 0.;
280 }
281
282 //________________________________________________________
283 void AliStrLine::InitTwoPoints(Double_t *pA, Double_t *pB){
284   // Initialization from the coordinates of two
285   // points in the space
286   Double_t cd[3];
287   for(Int_t i=0;i<3;i++)cd[i] = pB[i]-pA[i];
288   InitDirection(pA,cd);
289 }
290
291 //________________________________________________________
292 AliStrLine::~AliStrLine() {
293   // destructor
294   if(fWMatrix)delete [] fWMatrix;
295 }
296
297 //________________________________________________________
298 void AliStrLine::PrintStatus() const {
299   // Print current status
300   cout <<"=======================================================\n";
301   cout <<"Direction cosines: ";
302   for(Int_t i=0;i<3;i++)cout <<fCd[i]<<"; ";
303   cout <<endl;
304   cout <<"Known point: ";
305   for(Int_t i=0;i<3;i++)cout <<fP0[i]<<"; ";
306   cout <<endl;
307   cout <<"Error on known point: ";
308   for(Int_t i=0;i<3;i++)cout <<TMath::Sqrt(fSigma2P0[i])<<"; ";
309   cout <<endl;
310   cout <<"Current value for the parameter: "<<fTpar<<endl;
311 }
312
313 //________________________________________________________
314 Int_t AliStrLine::IsParallelTo(AliStrLine *line) const {
315   // returns 1 if lines are parallel, 0 if not paralel
316   Double_t cd2[3];
317   line->GetCd(cd2);
318   Double_t vecpx=fCd[1]*cd2[2]-fCd[2]*cd2[1];
319   if(vecpx!=0) return 0;
320   Double_t vecpy=-fCd[0]*cd2[2]+fCd[2]*cd2[0];
321   if(vecpy!=0) return 0;
322   Double_t vecpz=fCd[0]*cd2[1]-fCd[1]*cd2[0];
323   if(vecpz!=0) return 0;
324   return 1;
325 }
326 //________________________________________________________
327 Int_t AliStrLine::Crossrphi(AliStrLine *line){
328   // Cross 2 lines in the X-Y plane
329   Double_t p2[3];
330   Double_t cd2[3];
331   line->GetP0(p2);
332   line->GetCd(cd2);
333   Double_t a=fCd[0];
334   Double_t b=-cd2[0];
335   Double_t c=p2[0]-fP0[0];
336   Double_t d=fCd[1];
337   Double_t e=-cd2[1];
338   Double_t f=p2[1]-fP0[1];
339   Double_t deno = a*e-b*d;
340   Int_t retcode = 0;
341   if(deno != 0.) {
342     fTpar = (c*e-b*f)/deno;
343   }
344   else {
345     retcode = -1;
346   }
347   return retcode;
348 }
349
350 //________________________________________________________
351 Int_t AliStrLine::CrossPoints(AliStrLine *line, Double_t *point1, Double_t *point2){
352   // Looks for the crossing point estimated starting from the
353   // DCA segment
354   Double_t p2[3];
355   Double_t cd2[3];
356   line->GetP0(p2);
357   line->GetCd(cd2);
358   Int_t i;
359   Double_t k1 = 0;
360   Double_t k2 = 0;
361   Double_t a11 = 0;
362   for(i=0;i<3;i++){
363     k1+=(fP0[i]-p2[i])*fCd[i];
364     k2+=(fP0[i]-p2[i])*cd2[i];
365     a11+=fCd[i]*cd2[i];
366   }
367   Double_t a22 = -a11;
368   Double_t a21 = 0;
369   Double_t a12 = 0;
370   for(i=0;i<3;i++){
371     a21+=cd2[i]*cd2[i];
372     a12-=fCd[i]*fCd[i];
373   }
374   Double_t deno = a11*a22-a21*a12;
375   if(deno == 0.) return -1;
376   fTpar = (a11*k2-a21*k1) / deno;
377   Double_t par2 = (k1*a22-k2*a12) / deno;
378   line->SetPar(par2);
379   GetCurrentPoint(point1);
380   line->GetCurrentPoint(point2);
381   return 0;
382 }
383 //________________________________________________________________
384 Int_t AliStrLine::Cross(AliStrLine *line, Double_t *point){
385
386   //Finds intersection between lines
387   Double_t point1[3];
388   Double_t point2[3];
389   Int_t retcod=CrossPoints(line,point1,point2);
390   if(retcod==0){
391     for(Int_t i=0;i<3;i++)point[i]=(point1[i]+point2[i])/2.;
392     return 0;
393   }else{
394     return -1;
395   }
396 }
397
398 //___________________________________________________________
399 Double_t AliStrLine::GetDCA(AliStrLine *line) const{
400   //Returns the distance of closest approach between two lines
401   Double_t p2[3];
402   Double_t cd2[3];
403   line->GetP0(p2);
404   line->GetCd(cd2);
405   Int_t i;
406   Int_t ispar=IsParallelTo(line);
407   if(ispar){
408     Double_t dist1q=0,dist2=0,mod=0;
409     for(i=0;i<3;i++){
410       dist1q+=(fP0[i]-p2[i])*(fP0[i]-p2[i]);
411       dist2+=(fP0[i]-p2[i])*fCd[i];
412       mod+=fCd[i]*fCd[i];
413     }
414     if(mod!=0){
415       dist2/=mod;
416       return TMath::Sqrt(dist1q-dist2*dist2);
417     }else{return -1;}
418   }else{
419      Double_t perp[3];
420      perp[0]=fCd[1]*cd2[2]-fCd[2]*cd2[1];
421      perp[1]=-fCd[0]*cd2[2]+fCd[2]*cd2[0];
422      perp[2]=fCd[0]*cd2[1]-fCd[1]*cd2[0];
423      Double_t mod=0,dist=0;
424      for(i=0;i<3;i++){
425        mod+=perp[i]*perp[i];
426        dist+=(fP0[i]-p2[i])*perp[i];
427      }
428      mod=sqrt(mod);
429      if(mod!=0){
430        dist/=mod;
431        return TMath::Abs(dist);
432      }else{return -1;}
433   }
434 }
435 //________________________________________________________
436 void AliStrLine::GetCurrentPoint(Double_t *point) const {
437   // Fills the array point with the current value on the line
438   for(Int_t i=0;i<3;i++)point[i]=fP0[i]+fCd[i]*fTpar;
439 }
440
441 //________________________________________________________
442 Double_t AliStrLine::GetDistFromPoint(Double_t *point) const {
443   // computes distance from point 
444   AliStrLine tmp(point,(Double_t *)fCd,kFALSE);
445   return this->GetDCA(&tmp);
446 }