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   for(Int_t i=0;i<3;i++){
 208     fP0[i]=source.fP0[i];
 209     fSigma2P0[i]=source.fSigma2P0[i];
 210     fCd[i]=source.fCd[i];
 211   }
 212   if(source.fWMatrix){
 213     fWMatrix = new Double_t [6];
 214     for(Int_t i=0;i<6;i++)fWMatrix[i]=source.fWMatrix[i];
 215   }
 216 }
 217
 218 //________________________________________________________
 219 AliStrLine& AliStrLine::operator=(const AliStrLine& source){
 220   // Assignment operator
 221   if(this !=&source){
 222     this->~AliStrLine();
 223     new(this)AliStrLine(source);
 224   }
 225   return *this;
 226 }
 227
 228 //________________________________________________________
 229 void AliStrLine::GetWMatrix(Double_t *wmat)const {
 230 // Getter for weighting matrix, as a [9] dim. array
 231   if(!fWMatrix)return;
 232   Int_t k = 0;
 233   for(Int_t i=0;i<3;i++){
 234     for(Int_t j=0;j<3;j++){
 235       if(j>=i){
 236         wmat[3*i+j]=fWMatrix[k++];
 237       }
 238       else{
 239         wmat[3*i+j]=wmat[3*j+i];
 240       }
 241     }
 242   }
 243 }
 244
 245 //________________________________________________________
 246 void AliStrLine::SetWMatrix(const Double_t *wmat) {
 247 // Setter for weighting matrix, strating from a [9] dim. array
 248   if(fWMatrix)delete [] fWMatrix;
 249   fWMatrix = new Double_t [6];
 250   Int_t k = 0;
 251   for(Int_t i=0;i<3;i++){
 252     for(Int_t j=0;j<3;j++)if(j>=i)fWMatrix[k++]=wmat[3*i+j];
 253   }
 254 }
 255
 256 //________________________________________________________
 257 void AliStrLine::InitDirection(Double_t *point, Double_t *cd){
 258   // Initialization from a point and a direction
 259   Double_t norm = 0.;
 260   for(Int_t i=0;i<3;i++)norm+=cd[i]*cd[i];
 261   if(norm) {
 262     norm = TMath::Sqrt(norm);
 263     for(Int_t i=0;i<3;i++) cd[i]/=norm;
 264   }
 265   else {
 266     Error("AliStrLine","Null direction cosines!!!");
 267   }
 268   SetP0(point);
 269   SetCd(cd);
 270   fTpar = 0.;
 271 }
 272
 273 //________________________________________________________
 274 void AliStrLine::InitTwoPoints(Double_t *pA, Double_t *pB){
 275   // Initialization from the coordinates of two
 276   // points in the space
 277   Double_t cd[3];
 278   for(Int_t i=0;i<3;i++)cd[i] = pB[i]-pA[i];
 279   InitDirection(pA,cd);
 280 }
 281
 282 //________________________________________________________
 283 AliStrLine::~AliStrLine() {
 284   // destructor
 285   if(fWMatrix)delete [] fWMatrix;
 286 }
 287
 288 //________________________________________________________
 289 void AliStrLine::PrintStatus() const {
 290   // Print current status
 291   cout <<"=======================================================\n";
 292   cout <<"Direction cosines: ";
 293   for(Int_t i=0;i<3;i++)cout <<fCd[i]<<"; ";
 294   cout <<endl;
 295   cout <<"Known point: ";
 296   for(Int_t i=0;i<3;i++)cout <<fP0[i]<<"; ";
 297   cout <<endl;
 298   cout <<"Error on known point: ";
 299   for(Int_t i=0;i<3;i++)cout <<TMath::Sqrt(fSigma2P0[i])<<"; ";
 300   cout <<endl;
 301   cout <<"Current value for the parameter: "<<fTpar<<endl;
 302 }
 303
 304 //________________________________________________________
 305 Int_t AliStrLine::IsParallelTo(AliStrLine *line) const {
 306   // returns 1 if lines are parallel, 0 if not paralel
 307   Double_t cd2[3];
 308   line->GetCd(cd2);
 309   Double_t vecpx=fCd[1]*cd2[2]-fCd[2]*cd2[1];
 310   if(vecpx!=0) return 0;
 311   Double_t vecpy=-fCd[0]*cd2[2]+fCd[2]*cd2[0];
 312   if(vecpy!=0) return 0;
 313   Double_t vecpz=fCd[0]*cd2[1]-fCd[1]*cd2[0];
 314   if(vecpz!=0) return 0;
 315   return 1;
 316 }
 317 //________________________________________________________
 318 Int_t AliStrLine::Crossrphi(AliStrLine *line){
 319   // Cross 2 lines in the X-Y plane
 320   Double_t p2[3];
 321   Double_t cd2[3];
 322   line->GetP0(p2);
 323   line->GetCd(cd2);
 324   Double_t a=fCd[0];
 325   Double_t b=-cd2[0];
 326   Double_t c=p2[0]-fP0[0];
 327   Double_t d=fCd[1];
 328   Double_t e=-cd2[1];
 329   Double_t f=p2[1]-fP0[1];
 330   Double_t deno = a*e-b*d;
 331   Int_t retcode = 0;
 332   if(deno != 0.) {
 333     fTpar = (c*e-b*f)/deno;
 334   }
 335   else {
 336     retcode = -1;
 337   }
 338   return retcode;
 339 }
 340
 341 //________________________________________________________
 342 Int_t AliStrLine::CrossPoints(AliStrLine *line, Double_t *point1, Double_t *point2){
 343   // Looks for the crossing point estimated starting from the
 344   // DCA segment
 345   Double_t p2[3];
 346   Double_t cd2[3];
 347   line->GetP0(p2);
 348   line->GetCd(cd2);
 349   Int_t i;
 350   Double_t k1 = 0;
 351   Double_t k2 = 0;
 352   Double_t a11 = 0;
 353   for(i=0;i<3;i++){
 354     k1+=(fP0[i]-p2[i])*fCd[i];
 355     k2+=(fP0[i]-p2[i])*cd2[i];
 356     a11+=fCd[i]*cd2[i];
 357   }
 358   Double_t a22 = -a11;
 359   Double_t a21 = 0;
 360   Double_t a12 = 0;
 361   for(i=0;i<3;i++){
 362     a21+=cd2[i]*cd2[i];
 363     a12-=fCd[i]*fCd[i];
 364   }
 365   Double_t deno = a11*a22-a21*a12;
 366   if(deno == 0.) return -1;
 367   fTpar = (a11*k2-a21*k1) / deno;
 368   Double_t par2 = (k1*a22-k2*a12) / deno;
 369   line->SetPar(par2);
 370   GetCurrentPoint(point1);
 371   line->GetCurrentPoint(point2);
 372   return 0;
 373 }
 374 //________________________________________________________________
 375 Int_t AliStrLine::Cross(AliStrLine *line, Double_t *point){
 376
 377   //Finds intersection between lines
 378   Double_t point1[3];
 379   Double_t point2[3];
 380   Int_t retcod=CrossPoints(line,point1,point2);
 381   if(retcod==0){
 382     for(Int_t i=0;i<3;i++)point[i]=(point1[i]+point2[i])/2.;
 383     return 0;
 384   }else{
 385     return -1;
 386   }
 387 }
 388
 389 //___________________________________________________________
 390 Double_t AliStrLine::GetDCA(AliStrLine *line) const{
 391   //Returns the distance of closest approach between two lines
 392   Double_t p2[3];
 393   Double_t cd2[3];
 394   line->GetP0(p2);
 395   line->GetCd(cd2);
 396   Int_t i;
 397   Int_t ispar=IsParallelTo(line);
 398   if(ispar){
 399     Double_t dist1q=0,dist2=0,mod=0;
 400     for(i=0;i<3;i++){
 401       dist1q+=(fP0[i]-p2[i])*(fP0[i]-p2[i]);
 402       dist2+=(fP0[i]-p2[i])*fCd[i];
 403       mod+=fCd[i]*fCd[i];
 404     }
 405     if(mod!=0){
 406       dist2/=mod;
 407       return TMath::Sqrt(dist1q-dist2*dist2);
 408     }else{return -1;}
 409   }else{
 410      Double_t perp[3];
 411      perp[0]=fCd[1]*cd2[2]-fCd[2]*cd2[1];
 412      perp[1]=-fCd[0]*cd2[2]+fCd[2]*cd2[0];
 413      perp[2]=fCd[0]*cd2[1]-fCd[1]*cd2[0];
 414      Double_t mod=0,dist=0;
 415      for(i=0;i<3;i++){
 416        mod+=perp[i]*perp[i];
 417        dist+=(fP0[i]-p2[i])*perp[i];
 418      }
 419      mod=sqrt(mod);
 420      if(mod!=0){
 421        dist/=mod;
 422        return TMath::Abs(dist);
 423      }else{return -1;}
 424   }
 425 }
 426 //________________________________________________________
 427 void AliStrLine::GetCurrentPoint(Double_t *point) const {
 428   // Fills the array point with the current value on the line
 429   for(Int_t i=0;i<3;i++)point[i]=fP0[i]+fCd[i]*fTpar;
 430 }
 431
 432 //________________________________________________________
 433 Double_t AliStrLine::GetDistFromPoint(Double_t *point) const {
 434   // computes distance from point
 435   AliStrLine tmp(point,(Double_t *)fCd,kFALSE);
 436   return this->GetDCA(&tmp);
 437 }