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   const Double_t prec=1e-14;
 317   Double_t cd2[3];
 318   line->GetCd(cd2);
 319   Double_t vecpx=fCd[1]*cd2[2]-fCd[2]*cd2[1];
 320   if(TMath::Abs(vecpx) > prec) return 0;
 321   Double_t vecpy=-fCd[0]*cd2[2]+fCd[2]*cd2[0];
 322   if(TMath::Abs(vecpy) > prec) return 0;
 323   Double_t vecpz=fCd[0]*cd2[1]-fCd[1]*cd2[0];
 324   if(TMath::Abs(vecpz) > prec) return 0;
 325   return 1;
 326 }
 327 //________________________________________________________
 328 Int_t AliStrLine::Crossrphi(AliStrLine *line){
 329   // Cross 2 lines in the X-Y plane
 330   const Double_t prec=1e-14;
 331   const Double_t big=1e20;
 332   Double_t p2[3];
 333   Double_t cd2[3];
 334   line->GetP0(p2);
 335   line->GetCd(cd2);
 336   Double_t a=fCd[0];
 337   Double_t b=-cd2[0];
 338   Double_t c=p2[0]-fP0[0];
 339   Double_t d=fCd[1];
 340   Double_t e=-cd2[1];
 341   Double_t f=p2[1]-fP0[1];
 342   Double_t deno = a*e-b*d;
 343   Int_t retcode = 0;
 344   if(TMath::Abs(deno) > prec) {
 345     fTpar = (c*e-b*f)/deno;
 346   }
 347   else {
 348     fTpar = big;
 349     retcode = -1;
 350   }
 351   return retcode;
 352 }
 353
 354 //________________________________________________________
 355 Int_t AliStrLine::CrossPoints(AliStrLine *line, Double_t *point1, Double_t *point2){
 356   // Looks for the crossing point estimated starting from the
 357   // DCA segment
 358   const Double_t prec=1e-14;
 359   Double_t p2[3];
 360   Double_t cd2[3];
 361   line->GetP0(p2);
 362   line->GetCd(cd2);
 363   Int_t i;
 364   Double_t k1 = 0;
 365   Double_t k2 = 0;
 366   Double_t a11 = 0;
 367   for(i=0;i<3;i++){
 368     k1+=(fP0[i]-p2[i])*fCd[i];
 369     k2+=(fP0[i]-p2[i])*cd2[i];
 370     a11+=fCd[i]*cd2[i];
 371   }
 372   Double_t a22 = -a11;
 373   Double_t a21 = 0;
 374   Double_t a12 = 0;
 375   for(i=0;i<3;i++){
 376     a21+=cd2[i]*cd2[i];
 377     a12-=fCd[i]*fCd[i];
 378   }
 379   Double_t deno = a11*a22-a21*a12;
 380   if(TMath::Abs(deno) < prec) return -1;
 381   fTpar = (a11*k2-a21*k1) / deno;
 382   Double_t par2 = (k1*a22-k2*a12) / deno;
 383   line->SetPar(par2);
 384   GetCurrentPoint(point1);
 385   line->GetCurrentPoint(point2);
 386   return 0;
 387 }
 388 //________________________________________________________________
 389 Int_t AliStrLine::Cross(AliStrLine *line, Double_t *point){
 390
 391   //Finds intersection between lines
 392   Double_t point1[3];
 393   Double_t point2[3];
 394   Int_t retcod=CrossPoints(line,point1,point2);
 395   if(retcod==0){
 396     for(Int_t i=0;i<3;i++)point[i]=(point1[i]+point2[i])/2.;
 397     return 0;
 398   }else{
 399     return -1;
 400   }
 401 }
 402
 403 //___________________________________________________________
 404 Double_t AliStrLine::GetDCA(AliStrLine *line) const{
 405   //Returns the distance of closest approach between two lines
 406   const Double_t prec=1e-14;
 407   Double_t p2[3];
 408   Double_t cd2[3];
 409   line->GetP0(p2);
 410   line->GetCd(cd2);
 411   Int_t i;
 412   Int_t ispar=IsParallelTo(line);
 413   if(ispar){
 414     Double_t dist1q=0,dist2=0,mod=0;
 415     for(i=0;i<3;i++){
 416       dist1q+=(fP0[i]-p2[i])*(fP0[i]-p2[i]);
 417       dist2+=(fP0[i]-p2[i])*fCd[i];
 418       mod+=fCd[i]*fCd[i];
 419     }
 420     if(TMath::Abs(mod) > prec){
 421       dist2/=mod;
 422       return TMath::Sqrt(dist1q-dist2*dist2);
 423     }else{return -1;}
 424   }else{
 425      Double_t perp[3];
 426      perp[0]=fCd[1]*cd2[2]-fCd[2]*cd2[1];
 427      perp[1]=-fCd[0]*cd2[2]+fCd[2]*cd2[0];
 428      perp[2]=fCd[0]*cd2[1]-fCd[1]*cd2[0];
 429      Double_t mod=0,dist=0;
 430      for(i=0;i<3;i++){
 431        mod+=perp[i]*perp[i];
 432        dist+=(fP0[i]-p2[i])*perp[i];
 433      }
 434      mod=sqrt(mod);
 435      if(TMath::Abs(mod) > prec){
 436        dist/=mod;
 437        return TMath::Abs(dist);
 438      }else{return -1;}
 439   }
 440 }
 441 //________________________________________________________
 442 void AliStrLine::GetCurrentPoint(Double_t *point) const {
 443   // Fills the array point with the current value on the line
 444   for(Int_t i=0;i<3;i++)point[i]=fP0[i]+fCd[i]*fTpar;
 445 }
 446
 447 //________________________________________________________
 448 Double_t AliStrLine::GetDistFromPoint(Double_t *point) const {
 449   // computes distance from point
 450   AliStrLine tmp(point,(Double_t *)fCd,kFALSE);
 451   return this->GetDCA(&tmp);
 452 }