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   fTpar(0),
  37   fDebug(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   fTpar(0),
  51   fDebug(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++) fSigma2P0[i] = 0.;
  60   if(twopoints){
  61     InitTwoPoints(point,cd);
  62   }
  63   else {
  64     InitDirection(point,cd);
  65   }
  66 }
  67
  68 //________________________________________________________
  69 AliStrLine::AliStrLine(Float_t *pointf, Float_t *cdf,Bool_t twopoints) :
  70   TObject(),
  71   fTpar(0),
  72   fDebug(0)
  73 {
  74   // Standard constructor - with float arguments
  75   // if twopoints is true:  point and cd are the 3D coordinates of
  76   //                        two points defininig the straight line
  77   // if twopoint is false: point represents the 3D coordinates of a point
  78   //                       belonging to the straight line and cd is the
  79   //                       direction in space
  80   Double_t point[3];
  81   Double_t cd[3];
  82   for(Int_t i=0;i<3;i++){
  83     point[i] = pointf[i];
  84     cd[i] = cdf[i];
  85     fSigma2P0[i] = 0.;
  86   }
  87   if(twopoints){
  88     InitTwoPoints(point,cd);
  89   }
  90   else {
  91     InitDirection(point,cd);
  92   }
  93 }
  94
  95 //________________________________________________________
  96 AliStrLine::AliStrLine(Double_t *point, Double_t *sig2point, Double_t *cd,Bool_t twopoints) :
  97   TObject(),
  98   fTpar(0),
  99   fDebug(0)
 100 {
 101   // Standard constructor
 102   // if twopoints is true:  point and cd are the 3D coordinates of
 103   //                        two points defininig the straight line
 104   // if twopoint is false: point represents the 3D coordinates of a point
 105   //                       belonging to the straight line and cd is the
 106   //                       direction in space
 107   for(Int_t i=0;i<3;i++) fSigma2P0[i] = sig2point[i];
 108   if(twopoints){
 109     InitTwoPoints(point,cd);
 110   }
 111   else {
 112     InitDirection(point,cd);
 113   }
 114 }
 115
 116 //________________________________________________________
 117 AliStrLine::AliStrLine(Float_t *pointf, Float_t *sig2point, Float_t *cdf,Bool_t twopoints) :
 118   TObject(),
 119   fTpar(0),
 120   fDebug(0)
 121 {
 122   // Standard constructor - with float arguments
 123   // if twopoints is true:  point and cd are the 3D coordinates of
 124   //                        two points defininig the straight line
 125   // if twopoint is false: point represents the 3D coordinates of a point
 126   //                       belonging to the straight line and cd is the
 127   //                       direction in space
 128   Double_t point[3];
 129   Double_t cd[3];
 130   for(Int_t i=0;i<3;i++){
 131     point[i] = pointf[i];
 132     cd[i] = cdf[i];
 133     fSigma2P0[i] = sig2point[i];
 134   }
 135   if(twopoints){
 136     InitTwoPoints(point,cd);
 137   }
 138   else {
 139     InitDirection(point,cd);
 140   }
 141 }
 142 //________________________________________________________
 143 void AliStrLine::InitDirection(Double_t *point, Double_t *cd){
 144   // Initialization from a point and a direction
 145   Double_t norm = 0.;
 146   for(Int_t i=0;i<3;i++)norm+=cd[i]*cd[i];
 147   if(norm) {
 148     norm = TMath::Sqrt(norm);
 149     for(Int_t i=0;i<3;i++) cd[i]/=norm;
 150   }
 151   else {
 152     Error("AliStrLine","Null direction cosines!!!");
 153   }
 154   SetP0(point);
 155   SetCd(cd);
 156   fTpar = 0.;
 157   SetDebug();
 158 }
 159
 160 //________________________________________________________
 161 void AliStrLine::InitTwoPoints(Double_t *pA, Double_t *pB){
 162   // Initialization from the coordinates of two
 163   // points in the space
 164   Double_t cd[3];
 165   for(Int_t i=0;i<3;i++)cd[i] = pB[i]-pA[i];
 166   InitDirection(pA,cd);
 167 }
 168
 169 //________________________________________________________
 170 AliStrLine::~AliStrLine() {
 171   // destructor
 172 }
 173
 174 //________________________________________________________
 175 void AliStrLine::PrintStatus() const {
 176   // Print current status
 177   cout <<"=======================================================\n";
 178   cout <<"Direction cosines: ";
 179   for(Int_t i=0;i<3;i++)cout <<fCd[i]<<"; ";
 180   cout <<endl;
 181   cout <<"Known point: ";
 182   for(Int_t i=0;i<3;i++)cout <<fP0[i]<<"; ";
 183   cout <<endl;
 184   cout <<"Error on known point: ";
 185   for(Int_t i=0;i<3;i++)cout <<TMath::Sqrt(fSigma2P0[i])<<"; ";
 186   cout <<endl;
 187   cout <<"Current value for the parameter: "<<fTpar<<endl;
 188   cout <<" Debug flag: "<<fDebug<<endl;
 189 }
 190
 191 //________________________________________________________
 192 Int_t AliStrLine::IsParallelTo(AliStrLine *line) const {
 193   // returns 1 if lines are parallel, 0 if not paralel
 194   Double_t cd2[3];
 195   line->GetCd(cd2);
 196   Double_t vecpx=fCd[1]*cd2[2]-fCd[2]*cd2[1];
 197   if(vecpx!=0) return 0;
 198   Double_t vecpy=-fCd[0]*cd2[2]+fCd[2]*cd2[0];
 199   if(vecpy!=0) return 0;
 200   Double_t vecpz=fCd[0]*cd2[1]-fCd[1]*cd2[0];
 201   if(vecpz!=0) return 0;
 202   return 1;
 203 }
 204 //________________________________________________________
 205 Int_t AliStrLine::Crossrphi(AliStrLine *line){
 206   // Cross 2 lines in the X-Y plane
 207   Double_t p2[3];
 208   Double_t cd2[3];
 209   line->GetP0(p2);
 210   line->GetCd(cd2);
 211   Double_t a=fCd[0];
 212   Double_t b=-cd2[0];
 213   Double_t c=p2[0]-fP0[0];
 214   Double_t d=fCd[1];
 215   Double_t e=-cd2[1];
 216   Double_t f=p2[1]-fP0[1];
 217   Double_t deno = a*e-b*d;
 218   Int_t retcode = 0;
 219   if(deno != 0.) {
 220     fTpar = (c*e-b*f)/deno;
 221   }
 222   else {
 223     retcode = -1;
 224   }
 225   return retcode;
 226 }
 227
 228 //________________________________________________________
 229 Int_t AliStrLine::CrossPoints(AliStrLine *line, Double_t *point1, Double_t *point2){
 230   // Looks for the crossing point estimated starting from the
 231   // DCA segment
 232   Double_t p2[3];
 233   Double_t cd2[3];
 234   line->GetP0(p2);
 235   line->GetCd(cd2);
 236   Int_t i;
 237   Double_t k1 = 0;
 238   for(i=0;i<3;i++)k1+=(fP0[i]-p2[i])*fCd[i];
 239   Double_t k2 = 0;
 240   for(i=0;i<3;i++)k2+=(fP0[i]-p2[i])*cd2[i];
 241   Double_t a11 = 0;
 242   for(i=0;i<3;i++)a11+=fCd[i]*cd2[i];
 243   Double_t a22 = -a11;
 244   Double_t a21 = 0;
 245   for(i=0;i<3;i++)a21+=cd2[i]*cd2[i];
 246   Double_t a12 = 0;
 247   for(i=0;i<3;i++)a12-=fCd[i]*fCd[i];
 248   Double_t deno = a11*a22-a21*a12;
 249   if(deno == 0.) return -1;
 250   fTpar = (a11*k2-a21*k1) / deno;
 251   Double_t par2 = (k1*a22-k2*a12) / deno;
 252   line->SetPar(par2);
 253   GetCurrentPoint(point1);
 254   line->GetCurrentPoint(point2);
 255   return 0;
 256 }
 257 //________________________________________________________________
 258 Int_t AliStrLine::Cross(AliStrLine *line, Double_t *point){
 259
 260   //Finds intersection between lines
 261   Double_t point1[3];
 262   Double_t point2[3];
 263   Int_t retcod=CrossPoints(line,point1,point2);
 264   if(retcod==0){
 265     for(Int_t i=0;i<3;i++)point[i]=(point1[i]+point2[i])/2.;
 266     return 0;
 267   }else{
 268     return -1;
 269   }
 270 }
 271
 272 //___________________________________________________________
 273 Double_t AliStrLine::GetDCA(AliStrLine *line) const{
 274   //Returns the distance of closest approach between two lines
 275   Double_t p2[3];
 276   Double_t cd2[3];
 277   line->GetP0(p2);
 278   line->GetCd(cd2);
 279   Int_t i;
 280   Int_t ispar=IsParallelTo(line);
 281   if(ispar){
 282     Double_t dist1q=0,dist2=0,mod=0;
 283     for(i=0;i<3;i++){
 284       dist1q+=(fP0[i]-p2[i])*(fP0[i]-p2[i]);
 285       dist2+=(fP0[i]-p2[i])*fCd[i];
 286       mod+=fCd[i]*fCd[i];
 287     }
 288     if(mod!=0){
 289       dist2/=mod;
 290       return TMath::Sqrt(dist1q-dist2*dist2);
 291     }else{return -1;}
 292   }else{
 293      Double_t perp[3];
 294      perp[0]=fCd[1]*cd2[2]-fCd[2]*cd2[1];
 295      perp[1]=-fCd[0]*cd2[2]+fCd[2]*cd2[0];
 296      perp[2]=fCd[0]*cd2[1]-fCd[1]*cd2[0];
 297      Double_t mod=0,dist=0;
 298      for(i=0;i<3;i++){
 299        mod+=perp[i]*perp[i];
 300        dist+=(fP0[i]-p2[i])*perp[i];
 301      }
 302      mod=sqrt(mod);
 303      if(mod!=0){
 304        dist/=mod;
 305        return TMath::Abs(dist);
 306      }else{return -1;}
 307   }
 308 }
 309 //________________________________________________________
 310 void AliStrLine::GetCurrentPoint(Double_t *point) const {
 311   // Fills the array point with the current value on the line
 312   for(Int_t i=0;i<3;i++)point[i]=fP0[i]+fCd[i]*fTpar;
 313 }
 314
 315 //________________________________________________________
 316 Double_t AliStrLine::GetDistFromPoint(Double_t *point) const {
 317   // computes distance from point
 318   AliStrLine tmp(point,(Double_t *)fCd,kFALSE);
 319   return this->GetDCA(&tmp);
 320 }