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