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