STEER/ESD/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 <TMath.h>
  27
  28 #include "AliStrLine.h"
  29
  30 using std::endl;
  31 using std::cout;
  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(const Double_t *const point, const Double_t *const 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   else
  67     InitDirection(point,cd);
  68
  69   SetIdPoints(id1,id2);
  70 }
  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   TObject(),
  76   fWMatrix(0),
  77   fTpar(0)
  78 {
  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   for(Int_t i=0;i<3;i++)
  86     fSigma2P0[i] = sig2point[i];
  87
  88   if(twopoints)
  89     InitTwoPoints(point,cd);
  90   else
  91     InitDirection(point,cd);
  92
  93   SetIdPoints(id1,id2);
  94 }
  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   TObject(),
  99   fWMatrix(0),
 100   fTpar(0)
 101 {
 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   fWMatrix = new Double_t [6];
 110   for(Int_t i=0;i<3;i++){
 111     fSigma2P0[i] = sig2point[i];
 112     for(Int_t j=0;j<3;j++)if(j>=i)fWMatrix[k++]=wmat[3*i+j];
 113   }
 114   if(twopoints)
 115     InitTwoPoints(point,cd);
 116   else
 117     InitDirection(point,cd);
 118
 119   SetIdPoints(id1,id2);
 120 }
 121
 122 //________________________________________________________
 123 AliStrLine::AliStrLine(const AliStrLine &source):
 124   TObject(source),
 125   fWMatrix(0),
 126   fTpar(source.fTpar)
 127 {
 128   //
 129   // copy constructor
 130   //
 131   for(Int_t i=0;i<3;i++){
 132     fP0[i]=source.fP0[i];
 133     fSigma2P0[i]=source.fSigma2P0[i];
 134     fCd[i]=source.fCd[i];
 135   }
 136   if(source.fWMatrix){
 137     fWMatrix = new Double_t [6];
 138     for(Int_t i=0;i<6;i++)fWMatrix[i]=source.fWMatrix[i];
 139   }
 140   for(Int_t i=0;i<2;i++) fIdPoint[i]=source.fIdPoint[i];
 141 }
 142
 143 //________________________________________________________
 144 AliStrLine& AliStrLine::operator=(const AliStrLine& source)
 145 {
 146   // Assignment operator
 147   if(this !=&source){
 148     TObject::operator=(source);
 149     for(Int_t i=0;i<3;i++){
 150       fP0[i]=source.fP0[i];
 151       fSigma2P0[i]=source.fSigma2P0[i];
 152       fCd[i]=source.fCd[i];
 153     }
 154
 155     delete [] fWMatrix;
 156     fWMatrix=0;
 157     if(source.fWMatrix){
 158       fWMatrix = new Double_t [6];
 159       for(Int_t i=0;i<6;i++)fWMatrix[i]=source.fWMatrix[i];
 160     }
 161     for(Int_t i=0;i<2;i++) fIdPoint[i]=source.fIdPoint[i];
 162   }
 163   return *this;
 164 }
 165
 166 //________________________________________________________
 167 void AliStrLine::GetWMatrix(Double_t *wmat)const {
 168 // Getter for weighting matrix, as a [9] dim. array
 169   if(!fWMatrix)return;
 170   Int_t k = 0;
 171   for(Int_t i=0;i<3;i++){
 172     for(Int_t j=0;j<3;j++){
 173       if(j>=i){
 174         wmat[3*i+j]=fWMatrix[k++];
 175       }
 176       else{
 177         wmat[3*i+j]=wmat[3*j+i];
 178       }
 179     }
 180   }
 181 }
 182
 183 //________________________________________________________
 184 void AliStrLine::SetWMatrix(const Double_t *wmat) {
 185 // Setter for weighting matrix, strating from a [9] dim. array
 186   if(fWMatrix)delete [] fWMatrix;
 187   fWMatrix = new Double_t [6];
 188   Int_t k = 0;
 189   for(Int_t i=0;i<3;i++){
 190     for(Int_t j=0;j<3;j++)if(j>=i)fWMatrix[k++]=wmat[3*i+j];
 191   }
 192 }
 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   Double_t norm = cd[0]*cd[0]+cd[1]*cd[1]+cd[2]*cd[2];
 199
 200   if(norm) {
 201     norm = TMath::Sqrt(1./norm);
 202     for(Int_t i=0;i<3;++i) {
 203       fP0[i]=point[i];
 204       fCd[i]=cd[i]*norm;
 205     }
 206     fTpar = 0.;
 207   }
 208   else AliFatal("Null direction cosines!!!");
 209 }
 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   Double_t cd[3];
 217   for(Int_t i=0;i<3;i++)cd[i] = pB[i]-pA[i];
 218   InitDirection(pA,cd);
 219 }
 220
 221 //________________________________________________________
 222 AliStrLine::~AliStrLine() {
 223   // destructor
 224   if(fWMatrix)delete [] fWMatrix;
 225 }
 226
 227 //________________________________________________________
 228 void AliStrLine::PrintStatus() const {
 229   // Print current status
 230   cout <<"=======================================================\n";
 231   cout <<"Direction cosines: ";
 232   for(Int_t i=0;i<3;i++)cout <<fCd[i]<<"; ";
 233   cout <<endl;
 234   cout <<"Known point: ";
 235   for(Int_t i=0;i<3;i++)cout <<fP0[i]<<"; ";
 236   cout <<endl;
 237   cout <<"Error on known point: ";
 238   for(Int_t i=0;i<3;i++)cout <<TMath::Sqrt(fSigma2P0[i])<<"; ";
 239   cout <<endl;
 240   cout <<"Current value for the parameter: "<<fTpar<<endl;
 241 }
 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   Double_t cd2[3];
 248   line->GetCd(cd2);
 249
 250   Double_t vecpx=fCd[1]*cd2[2]-fCd[2]*cd2[1];
 251   Double_t mod=TMath::Abs(fCd[1]*cd2[2])+TMath::Abs(fCd[2]*cd2[1]);
 252   if(TMath::Abs(vecpx) > prec*mod) return 0;
 253
 254   Double_t vecpy=-fCd[0]*cd2[2]+fCd[2]*cd2[0];
 255   mod=TMath::Abs(fCd[0]*cd2[2])+TMath::Abs(fCd[2]*cd2[0]);
 256   if(TMath::Abs(vecpy) > prec*mod) return 0;
 257
 258   Double_t vecpz=fCd[0]*cd2[1]-fCd[1]*cd2[0];
 259   mod=TMath::Abs(fCd[0]*cd2[1])+TMath::Abs(fCd[1]*cd2[0]);
 260   if(TMath::Abs(vecpz) > prec) return 0;
 261
 262   return 1;
 263 }
 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   Double_t p2[3];
 271   Double_t cd2[3];
 272   line->GetP0(p2);
 273   line->GetCd(cd2);
 274   Double_t a=fCd[0];
 275   Double_t b=-cd2[0];
 276   Double_t c=p2[0]-fP0[0];
 277   Double_t d=fCd[1];
 278   Double_t e=-cd2[1];
 279   Double_t f=p2[1]-fP0[1];
 280   Double_t deno = a*e-b*d;
 281   Double_t mod=TMath::Abs(a*e)+TMath::Abs(b*d);
 282   Int_t retcode = 0;
 283   if(TMath::Abs(deno) > prec*mod) {
 284     fTpar = (c*e-b*f)/deno;
 285   }
 286   else {
 287     fTpar = big;
 288     retcode = -1;
 289   }
 290   return retcode;
 291 }
 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   Double_t p2[3];
 299   Double_t cd2[3];
 300   line->GetP0(p2);
 301   line->GetCd(cd2);
 302   Int_t i;
 303   Double_t k1 = 0;
 304   Double_t k2 = 0;
 305   Double_t a11 = 0;
 306   for(i=0;i<3;i++){
 307     k1+=(fP0[i]-p2[i])*fCd[i];
 308     k2+=(fP0[i]-p2[i])*cd2[i];
 309     a11+=fCd[i]*cd2[i];
 310   }
 311   Double_t a22 = -a11;
 312   Double_t a21 = 0;
 313   Double_t a12 = 0;
 314   for(i=0;i<3;i++){
 315     a21+=cd2[i]*cd2[i];
 316     a12-=fCd[i]*fCd[i];
 317   }
 318   Double_t deno = a11*a22-a21*a12;
 319   Double_t mod = TMath::Abs(a11*a22)+TMath::Abs(a21*a12);
 320   if(TMath::Abs(deno) < prec*mod) return -1;
 321   fTpar = (a11*k2-a21*k1) / deno;
 322   Double_t par2 = (k1*a22-k2*a12) / deno;
 323   line->SetPar(par2);
 324   GetCurrentPoint(point1);
 325   line->GetCurrentPoint(point2);
 326   return 0;
 327 }
 328 //________________________________________________________________
 329 Int_t AliStrLine::Cross(AliStrLine *line, Double_t *point)
 330 {
 331
 332   //Finds intersection between lines
 333   Double_t point1[3];
 334   Double_t point2[3];
 335   Int_t retcod=CrossPoints(line,point1,point2);
 336   if(retcod==0){
 337     for(Int_t i=0;i<3;i++)point[i]=(point1[i]+point2[i])/2.;
 338     return 0;
 339   }else{
 340     return -1;
 341   }
 342 }
 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   Double_t p2[3];
 350   Double_t cd2[3];
 351   line->GetP0(p2);
 352   line->GetCd(cd2);
 353   Int_t i;
 354   Int_t ispar=IsParallelTo(line);
 355   if(ispar){
 356     Double_t dist1q=0,dist2=0,mod=0;
 357     for(i=0;i<3;i++){
 358       dist1q+=(fP0[i]-p2[i])*(fP0[i]-p2[i]);
 359       dist2+=(fP0[i]-p2[i])*fCd[i];
 360       mod+=fCd[i]*fCd[i];
 361     }
 362     if(TMath::Abs(mod) > prec){
 363       dist2/=mod;
 364       return TMath::Sqrt(dist1q-dist2*dist2);
 365     }else{return -1;}
 366   }else{
 367      Double_t perp[3];
 368      perp[0]=fCd[1]*cd2[2]-fCd[2]*cd2[1];
 369      perp[1]=-fCd[0]*cd2[2]+fCd[2]*cd2[0];
 370      perp[2]=fCd[0]*cd2[1]-fCd[1]*cd2[0];
 371      Double_t mod=0,dist=0;
 372      for(i=0;i<3;i++){
 373        mod+=perp[i]*perp[i];
 374        dist+=(fP0[i]-p2[i])*perp[i];
 375      }
 376      if(TMath::Abs(mod) > prec) {
 377        return TMath::Abs(dist/TMath::Sqrt(mod));
 378      } else return -1;
 379   }
 380 }
 381 //________________________________________________________
 382 void AliStrLine::GetCurrentPoint(Double_t *point) const {
 383   // Fills the array point with the current value on the line
 384   for(Int_t i=0;i<3;i++)point[i]=fP0[i]+fCd[i]*fTpar;
 385 }
 386
 387 //________________________________________________________
 388 Double_t AliStrLine::GetDistFromPoint(const Double_t *point) const
 389 {
 390   // computes distance from point
 391   AliStrLine tmpline(point, fCd, kFALSE);
 392   return GetDCA(&tmpline);
 393 }
 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   Double_t p1= fCd[0]*fP0[0]+fCd[1]*fP0[1];
 404   Double_t p2=fCd[0]*fCd[0]+fCd[1]*fCd[1];
 405   Double_t delta=p1*p1-p2*(fP0[0]*fP0[0]+fP0[1]*fP0[1]-r*r);
 406   if(delta<0.){
 407     t1=-1000.;
 408     t2=t1;
 409     return kFALSE;
 410   }
 411   delta=TMath::Sqrt(delta);
 412   t1=(-p1-delta)/p2;
 413   t2=(-p1+delta)/p2;
 414
 415   if(t2<t1){
 416     // use delta as a temporary buffer
 417     delta=t1;
 418     t1=t2;
 419     t2=delta;
 420   }
 421   if(TMath::AreEqualAbs(t1,t2,1.e-9))t1=t2;
 422   return kTRUE;
 423 }